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I.  Introduction 

In  this  report  a  study  has  bean  made  of  ths  induced  currents  on 
a  two-wire  (parallel-cylinder)  transmission  line  due  to  incident  field 
excitation.  The  basic  assumptions  are  that  the  wires  are  perfectly 
conducting,  parallel.  Infinitely  long  and  that  their  radii  are  very 
small  compared  to  the  wavelength  of  the  excitation.  Thus,  the  Induced 
currents  along  the  wires  are  oriented  in  the  direction  of  the  cylinder 
axis  of  each  wire,  they  are  independent  of  the  asimuthal  angle  around 
each  wire,  and  hence  we  may  replace  each  wire  by  an  equivalent  filament. 
Boundary  conditions  are  applied  at  the  wire  surface  so  that  wire  thick¬ 
ness  is  taken  into  account.  Previous  work  related  to  this  problem 
included  Schelkunof f 1 s  Laplace  transform  solution  [1]  and  Marin's 
transient  solution  [2]  of  the  two-wire  problem. 

The  formulation  of  the  problem  is  effected  by  using  the  spectral 
concept.  The  excitation  field  has  been  resolved  into  its  Fourier 
(spatial)  components  and  the  Induced  current  has  been  expressed  in  terms 
of  a  superposition  integral.  First,  the  case  of  plane  wave  excitation 
has  been  treated.  This  led  to  a  straight-forward  solution  of  the  integral 
representation  and  a  simple  solution  for  the  induced  currents.  Using 
the  results,  the  variations  of  the  induced  currents  have  been  shown  as 
a  function  of  the  distance  between  the  wires  (while  the  ratio  of  the 
distance  between  the  wires  to  the  radius  of  the  wire  has  been  kept  fixed). 
Also,  the  variation  of  the  currents  has  been  shown  as  a  function  of  the 
radius  of  the  wire  (while  the  distance  between  the  wires  has  been  kept 
fixed) . 

Special  attention  has  been  paid  to  the  case  of  voltage  source 
excitation.  In  this  case,  by  deforming  the  contour  of  integration  in 
the  induced  current  formulation  and  applying  Cauchy's  theorem,  one 
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obtains  the  solution  In  terns  of  the  so-called  leaky  modes,  the  TEM  mode, 
and  a  continuous  spectrum.  The  leaky  modes  are  not  proper  solutions  of 
Maxwell's  equation  and,  in  contrast  to  the  discrete  eigenmodes,  which 
exist  in  regions  of  finite  extent  bounded  by  impermeable  walls,  do  not 
possess  orthogonality  and  completeness  properties,  and  therefore  must 
be  supplemented  by  a  continuous  spectrum  of  characteristic  modes  to  permit 
the  representation  of  an  arbitrary  function.  Nevertheless,  despite 
their  physically  unacceptable  behavior  in  the  entire  domain,  they  can 
be  employed  to  obtain  a  convergent  representation  of  the  field  solution 
in  certain  regions. 

The  preceding  modal  solution  may  be  obtained  via  a  straight-forward 
approach  using  the  general  network  parameters  of  the  method  of  moments  [3]. 
Solving  for  the  zeros  of  the  total  Impedance  of  the  system  leads  as 
expected  to  the  same  equation  that  determines  the  poles  of  the  integrand 
in  the  integral  solution  for  the  induced  currents. 

The  influence  of  the  electrical  dimensions  of  the  structure  on  the 
induced  currents  and  the  relative  importance  of  the  higher  modes  are  examined 
by  concrete  numerical  examples.  By  looking  at  the  results,  it  can  be 
readily  seen,  as  one  can  expect,  that  the  attenuation  constants  of  the 
higher  modes  decrease  as  the  distance  between  the  wires  is  increased. 

In  contrast,  the  attenuation  constants  increase  for  fixed  distance  bet¬ 
ween  the  wires  as  the  radius  is  decreased.  The  magnitudes  of  the  TEM 
and  the  higher  order  modal  currents  have  also  been  compared  at  a  fixed 
point  along  the  z  axis.  As  one  proceeds  along  the  z  axis,  the  higher 
order  modes  are  attenuated  exponentially,  while  the  TEM  remains  constant. 

These  studies  apply  particularly  to  problems  involving  high 
frequency  excitation,  such  that  the  wire  separation  is  comparable  to 
or  larger  than  a  wavelength,  in  which  case  the  usual  TEM  mode  analysis 
is  inadequate. 
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II.  Formulation  of  the  Problem 


Consider  the  problem  of  two  perfectly  conducting,  parallel,  thin, 
infinitely  long  wire?  with  radii  a  separated  by  a  distance  d.  The  *- 
directed  induced  currents  on  the  wires  are  due  to  an  Incident  field 
excitation  which  may  be  either  symmetric  or  antisymmetric  aa  shown 
in  Figure  1.  The  symmetric  excitation  (which  conslata  of  two  plane  waves) 
yields  a  magnetic  wall  it  y  ■  Similarly,  the  antisymmetric  excitation 
yields  an  electric  wall  at  y  ■  An  arbitrary  plane  wa#e  can  be 
represented  in  terms  of  a  sum  of  symmetric  and  antisymmetric  excitations. 

We  assume  that  j  «  1  and  y  «  1,  where  X  is  the  wavelength.  Thus,  the 
induced  current  is  along  the  cylinder  axis  of  each  wire,  and  is  independent 
of  the  azimuthal  angle  around  each  wire. 

The  magnetic  vector  potential  for  this  problem  is  given  by  (for 
exp(jd't)  time  convention) 

A  (x.y.z)  -  *'f(x,y,t)  (1) 


where 

't'+(x,y,z  ) 


•  | 


e-Jk/2  +  y2  +■  f.  -  e’)2 
4 it/x2  +  y2  +  (  z  -  t')* 


e-jkv42  +  (y  -  d)2  +  (z  -  z’)2 
+  (z  -  z’)2 


dz' 


(2) 


I(z)  denotes  the  current  along  the  wire,  x,y,z,  are  the  coordinates  of 
a  field  point  in  a  rectangular  coordinate  system,  z'  is  the  coordinate 
of  a  source  point,  k  is  the  wavenumber,  and  the  upper  and  lower  signs 
are  for  symmetric  and  antisymmetric  excitations,  respectively. 

The  Fourier  transform  of  i'(x,y,z)  with  respect  to  the  cylinder 
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axis  z  is. 


H'±(x,y,k2)  -  |  H±(x,y, 

—CO 


and  it  may  be  given  by. 


-jk  z 

z)e  dz 


(3) 


4±(x’y’V  =  Vf  W  tHo2>  (VP  ±  Ho2)<V2)]  <4> 

where  It(kz)  is  the  transform  of  I±(z);  kp  the  radial  wavenumber  is  defined 


as 


kp  -  (k2  -  k2>1/2 


(5) 


2  ^ 


2,1/2 


and  px  =  (x“  +  y")  ",  =  [x“  4  (y  -  d)"l 

<2) 

The  Hq  is  the  Hankel  function  of  the  second  kind  of  zero  order,  and 
we  have  used  the  identity  [4] 


-jk/p 


JF77 


1TT7 


2j 


(•n  ik7z 

Ho  (kpp)  e  dk-. 


«*> 


According  to  our  assumption  ^  <<  1.  Therefore  we  can  make  the  following 
approximation 

jr* 


(y  -  d)‘ 


2  2  2  :  d(J  "  d  +  i  £>>  :  d 

x  +  >‘  ■  a  '  d" 


(7) 


“?  '7  ‘7 

Using  (7)  )  (4)  way  be  simplified  for  x  +  y  —  n*"  ;uul  wo  obtain. 


V±(xfy,k2) 


!  "o2,V>' 


(«) 


The  z-directed  scattered  field  due  to  the  induced  currents  is  obtained 
from  V  according  to 

Ez±(x*y’2)  =  +  <x,y,Z>  (9) 

The  field  can  also  be  written  as  a  function  of  Y,  the  Kntirier  transform 
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of  4\  Using  the  inverse  transform  of  (3),  namely 


Y+(x,y,z)  -  —  'r  +  (x,y,k  )e  dk 


and  substituting  (10)  into  (9)  results  in 


E^(x,y,?.)  -  ^ 


7  k2  -  k2  jk  z 

|  7±(x,y,kz)e  z  dkz 


Hence,  the  transformed  solution  to  the  field  is  readily  given  by 
k2  -  k2 

¥z±(x*y*kz^  <1: 

where  V  is  given,  for  the  general  case,  by  (4). 

The  total  z-direeted  field,  with  the  conducting  cylinders  present, 
is  the  sum  of  the  incident  and  scattered  fields,  that  is, 

E  -  E1  +  ES 


For  each  cylinder  the  boundary  condition  E  *  0  must  be  satisfied  on  the 

z 

cylinder  surface.  Without  loss  of  generality,  we  can  deal  with  the  cylinder 
centered  at  the  origin.  Thus,  we  obtain 


E*  (x,y,z)  +  E®  (x,y,z) 


2  .  2  2  , 
x  4-  y  ■  a  ,  — 00  <  z  < 


The  transformed  solution  of  the  field  should  also  satisfy  (13),  that  is. 


E  (x,y,k  )  +  E  (x,y,k  ) 


z  z, 


2  2  2 
'x  +y  “a 


where 


E*  (x,y,k  ) 
z+  z 


(x.y, 

,  4 


-jk  z 

z)e  2  dz 


E  (x,y,k  )  is  given  by  (12). 
z  z 


Using  E  together  with  the  approximate  ¥ ,of  (8)  in  (14)  results  in  the 
z+  ± 


following  expression  for  the  transformed  current 

I  J 

M‘,>  ■  o»: 

Vo’V' 1  «“’<v*>i 

Applying  straight-forward  inverse  transfora  on  (16),  one  readily  obtains 
the  current  on  the  wire.  We  have, 

.  _2*  f' ,2  +  .2..2  jk_. 

i  k7  tt  j  2.  (2)  (2) - *  *  dk  07) 

V“o  <V> 1  "f  <v>i 

In  case  of  an  incident  wave  which  is  neither  symmetric  nor  antisymmetric, 

it  may  be  represented  as  a  superposition  of  symmetric  and  antisymmetric 

incident  waves  as  follows 

F*  -  E1  +  E1 
2  2+  2- 

and  the  currents  will  be  given  then  by 

Ix(z)  -  I+(z)  +  I_(z)  (lg) 

for  the  wire  centered  at  x  -  0,  y  -  0  and 

l2(z)  «  l+(z)  -  I_(Z)  (19) 

for  the  wire  centered  at  x  -  0,  y  -  d.  From  (17)  we  quote 


T+(*) 


I  3 


Ft.  (x’y*kz>  2  2  2 

•*'  lx  +  y  »  a 

^|Ho2>(kpa>  +  Ho2>(kpd>l 


e  2  dk 


I  <*)  - 


E,  <-->'-kz>  2  2  2 
_ .  ~ _ [x  t  y  »  a 

kX2><v>  - 


e  z  dk 


[II.  Plane  Wave  Excitation 

Figure  2  shows  a  plane  wave  incident  upon  the  two-wire  structure. 

The  electric  incident  field  has  been  divided  into  the  sum, 

Eln°  -  e5"C  +  EinC  (21) 

—  — h  — e 

The  subscripts  h  and  e  denote  H  and  E  polarisations,  respectively.  H 
modes  ate  those  for  which  E  is  parallel  to  the  plane  spanned  by  the  z  and  £. 
ft  represents  a  unit  vector  in  the  direction  of  propagation;  the  carat 
represents  a  unit  vector.  E  modes  are  those  for  which  E  is  normal  to  this 
plane.  Hence,  we  have 

Einc  .  [Einc  .  (k  x  z)]  _a_x.j_  .  (22) 

|K  x  $ | 2 

E*nC  -  Elnc  -  EinC  (23) 

— n  — e 

Since  we  have  already  limited  ourselves  to  a  thin  wire  problem  which  can 
be  replaced  with  a  filamentary  model,  the  E-polarized  component  does 
not  contribute  to  the  induced  current.  Therefore  we  will  consider  only 
the  H  component.  The  last  is  given  for  the  incident  wave  depicted  in 
Figure  2  by 

E^nC  (x,y,z)  -  E  (-x  cos  0  cos  -  y  cos  0  sin  $  + 

— h  o  o  o  o  o 

+1(k  +  k  y  +  k  z)  (24) 

~  ,  Q  .  xo  yo7  zo 

z  sin  R  )e 
o 

where  k  ,  k  ,  k  ,  the  axial  wave  numbers  are 
xo  yo’  zo 

k  *  k  sin  0  cos  d>  ,  k  *  sin  6  sin  <fc  ,  and  k  *  k  cos  0  ,  (25) 

xo  o  o  yo  o  o  zo  o 

and  0q,  <pQ  are,  respectively,  the  polar  and  azimuthal  angles  of  the 

incident  wave  in  spherical  coordinates. 

inc 

The  z  component  of  E  may  be  represented  as  a  sum  of  symetric 
and  antisymmetric  incident  waves  as  follows 


Einc  +  Einc 
2+ 


8 


fH 

h 

i.i 

4 » 


H 

r 

F 


Inc 

where  E  ,  the  symmetric  excitation,  ia  given  by 
♦ 

E^nc(x,y,t)  -  y  E  sin  0  [e^”*  +  ^  +  + 

i  o  o 

.1"W  +  kyo(<l  -  »>  +  k„*l 
inc 

and  ,  the  antisymmetric  excitation,  is  feiven  by 

,  j(k  x  +  k  y  +  k  z) 

E*nc(x,y,z)  ■  4  E  sin  0  |e  xo  yo  zo  - 

z  zoo 


ej[kxoX  +  kyo(d  -  y>  +  ‘‘zo2^ 


(27) 


(28) 


The  transformed  field  is  readily  obtained  by  applying  the  Fourier 
transform.  We  have 

tf*  -  B +  i8inc  (29) 


i 


II 

.  \ 
ti 

K 

I 


where 

_  jk  x  jk  y  jk  <d  -  y) 

Ez+(x,y,kz)  "  ,EoSln  0oe  X°  [e  y°  +  e  y°  l6(kz-  kJO)  (30) 

jk  x  jk  y  jV.  (d  -  y) 

Ez  (x,y,kz)  -  TTEosin  6^  *°  [e  y0  -  e  y0  ]6(kz-  k^)  (31) 

2  2  2  2 

where  the  6  denotes  theimpulse  function.  For  x  +  y  «  a«  A  , 


H 

K  - 

Li 


t 


t 


(30)  and  (31)  reduce,  respectively,  to. 


E  (x,y,k  ) 
=+  2 


and 


E  (x.y.k  ) 
z  z 


jk  d 

,  ,  ,  =  ttE  sin  0(1  +  e  yo  )  6(k  -  k  )  (32) 

x2  +  y2  -  a2  0  0  z  20 


jk  d 

2  2  2  ~  "V1"  6o(1  e  y°  >  6<k,  -  kz0)  <33> 

x  +  y  *  a 


substituting  (32)  and  (33)  in  (18),  (19)  and  (20)  yields 


2E 


V2>  -  * 


oO 


r  Jk  d 

_ 1  4  e  yo 

H(2^(k  a)  +  H(2) 
.  o  v  po  o 


+  — ; 


1  -  e 


jk  d 
yo 


q2) (k  d)  H(2)(k  a)  -  H(2)(k  d) 
o  po  o  po  o  po 


jk  z 
zo 


(34) 
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.i  I  ‘ 
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IV 


Solution  of  the  Superposition  Integral  in  Terms  of  the  TEM  Mod*  and  the 

Leaky  Modes 


The  Integral  representation  for  the  current  (Equation  (17J)  is  re¬ 
written  for  convenience  as 

QL> 

f+(k  ,x.y)  ,  9  e  *  dk  (37) 

L~  x"  +  y“  -  a  z 

—OP  •  v 


where 

f+  (k  *x,y)  -  j  2  2 

!x"  +  y  -  a 


e;  (X,y,k.) 

2  ,  2  2 

> 

x  +  y  •  a 

k2[H(2)(k  a)  t  H(2)(k  d)] 

C  O  CO  o  p 


The  singularities  of  the  integrand  are  two  branch  points  at  k^  ■  +k. 

On  Introducing  a  small  loss  perturbation  one  finds  that  for  positive 
values  In  the  z  direction  (z  >  0) ,  the  contour  runs  on  the  top  sheet 
of  the  kz  plane  along  the  real  axis,  below  the  branch  point  -k  and  above 
+k.  If  the  integrand  in  (37)  does  not  have  poles  in  the  upper  half 
of  the  top  sheet  of  the  k^  plane,  it  becomes  possible  to  deform  the 
original  path  of  integration,  the  real  axis,  from  to  +°°  ,  to  a 
new  path  ABCD  as  shown  in  Figure  3.  The  path  ABCD  may  be  modified  by 
letting  it  cross  the  branch  cut  twice,  entering  into  the  bottom  sheet 
in  this  process.  The  portion  of  the  path  in  the  lower  sheet,  namely 
C' ,  may  now  be  deformed  into  a  new  path  C",  as  shown  in  Figure  4,  with 
the  purpose  of  improving  the  convergence  of  the  integral  [2,  6,  7].  In 
the  process  of  carrying  out  the  last  step,  we  cross  the  poles  of  the 
Integrand  in  (37).  These  poles  are  solutions  of  the  equation 

H(2)  (/ i 2  -  k2a)  l  H(2)  (X2"-  k2d)  -  0  (38) 

o  z  o  z 

Equation  (38)  corresponds  to  Equation  (30)  of  Marin  (2].  Lying  in 
the  second  quadrant  of  the  k^  plane,  these  poles  are  associated  with 
fields  that  decay  as  they  propagate  in  the  z  direction.  HoVever,  since 
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-  TOP  SHEET 

- BOTTOM  SHEET 


Deforming  Che  contour  to  cross  the  leaky  poles. 


these  poles  btlong  to  ths  bottom  shoot  of  tho  k(  plono,  tho  fields 
associated  with  thorn  show  o  growth  bohovlor  in  tho  ttonovoroo  direction. 
For  thio  reeson,  thooo  nodes  oro  referred  to  oo  the  Improper,  or  leeky, 
modes  [7). 

In  thlo  fashion,  when  deforming  the  contour  no  contribution  will 
arise  from  eny  segment  of  the  Infinite  circle  (A  +  D)  on  the  proper  sheet 
of  the  complex  k(  plane.  On  spplying  Cauchy ' j  theorem,  one  hee 

I±(«)  ^  ^residues  of  Jf>  (f  (k^.x.y)  2  2  2 

at  the  leaky  poles  when  C'  la  deformed  to  C” .  )j 

f+(k  ,x,y)|  0  ,  2exp(jk  z)dk  (39) 

|  x“  +  y  ■  a 
C 

The  first  term  In  Equation  (39)  is  the  contribution  of  the  leaky 
poles,  and  it  may  be  rewritten  aa  (Appendix  Al) 


!.<*)  - 

1 leaky 


JL 

n  k, 


Hj,C(x,y.k_+) 


k _ *♦  *P  x-  4  y"  -  a"  «  p J Sp± * 


(40) 


where  k  ■  k  are  tho  consecutive  roots  of  the  equation 

*  *P  + 


«J2)(A2  -  k2  a)  +  hJ2><^?'-'  k2  d)  »  0 


and  will  be  referred  to  as  even  poles;  k^  *  k  are  the  eonsectuve 
roots  of  the  equation 


»<2)  (/F 


kr.  u) 


O  7.  '  0 

and  will  be  referred  to  as  odd  poles,  and 

k  -  (k2  -  k2  )1/2 
PP±  1  *pi' 


?  k2  d)  -  0 

7. 


In  the  above  equations  we  exclude,  of  course,  those  roots  for  which 
Re(k  )  <  -k  or  Re(k  )  >  0  and  Im  <k  <  0.  This  series.  In  (40\  is  often 
the  dominant  contribution  to  the  current  of  the  wires  [6\. 


for  negative  valuta  In  the  t  direction  (a  <  0) ,  tha  contour 
run*  on  tha  top  ahaat  of  tha  k#  plana  along  tha  real  axis,  above  tha 
branch  point  -k  and  balov  4k.  tf  tha  Integrand  In  (37)  doaa  not  hava 
polaa  in  tha  lover  half  of  tha  top  ahaat  of  the  k(  plane,  than  Mklng 
similar  contour  deformation  in  tha  lover  part  of  tha  k#  plana,  one  crosses 
tha  leaky  polaa,  vhich  are  tha  roota  of  Equation  (38)  that  lie  in  tha 
fourth  quadrant  and  are  the  image  of  tha  pravioua  polaa.  Hence  for 
t  <  0  vt  nay  apply  (AO)  replacing  k  by  -k  , , 

■pi  ipz 

The  aacond  term  ia  in  fact  tha  integral  along  a  branch  cut  etarting 
at  -k  and  parallel  to  tha  imaginary  axis  in  the  upper  half  of  the  k 

c 

plana  [8].  It  results  in  tha  contribution  of  the  branch  point  which, 
for  antisymmetric  excitation,  may  be  interpreted  as  the  TEM 
mode,  represented  by  (Appendix  A2) 


Ii«) 

TEM 


. _  -jrinc 

0  *n(£)  *- 


(x,y,-k) 


(Al) 


For  symmetric  excitation  the  branch  point  contribution  is  zero.  The 
other  contribution  is  the  branch  cut  contribution,  which  is  rather  time- 
consuming  to  evaluate,  but  is  negligible  [8]. 

Hence,  the  current  induced  in  the  case  of  an  arbitrary  incident 
wave  will  be  given  by 

I.(a>«I(z)  +  1. (a)  +  I  (z)  (A2) 

A  TEM  leaky  "leaky 

for  the  wire  centered  at  x  »  0,  y  ■  0,  and  by 

I  (z)  -  -I_  (z)  4  I  (*)  -  I_(z)  (A3) 

"TEM  leaky  "leaky 

for  the  wire  centered  at  x  ■  0,  y  »  d,  where  I.(z)  and  I  (z)  are  given 

leaky  "TEM 

by  Equation  (AO)  and  Equation  (Al),  respectively. 
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V.  Voltan®  Source  Excitation 


Consider  the  id®*?,  volt®*®  *®n®r®tor  c*nt®r®d  ®t  th®  ori*in  •• 
depleted  in  Figure  5.  Th®  equivalent  ®up®rpo®ition  of  aymetric  ®nd 
antisyommtric  excitations  i®  ®hovn  in  Fi*ur®  6.  For  thi®  c®®®,  th®  inci¬ 
dent  wave  i®  given  by  1 3] 

21 

j  vfi(«)  x"  +  y  -  r 

K*"'* 


slnc  .  fv5(,>  +  ,S  ‘  * 

1  I  o  ®ls®wh«r® 


(44) 


where  6  denotes  the  impulse  function.  EfnC  in  (44)  may  b®  caat  into  the 
form 


.inc 


„inc  _  a 

E.  < 


inc  +  pine 

v/2  6(«) 

x2  +  y2  •* 

v/2  6(a) 

x2  +  (y  - 

0 

1 

elsewhere 

v/2  6(a) 

x2  +  y2  ■ 

-v/2  6(a) 

x2  +  (y- 

i  o 

elsevhere 

(45) 

(46) 
2  2 


The  tranformed  field  may  be  readily  obtained  and  at  x  +  y 
we  have 

«!“<*.y.k,)  I  J  2  2  ■  <C  2  2  .  ,2 

*  x  +  y  -  a  +  17 


where 


(47) 

EinC(x.y.k  )  1  ,  ,  2  "  v/2 

*+  *  |x*  +  y  “  a 

(48) 

^nc(x,y.k  )  1  2  2  2  "  v/2 

*  z  |xz  +  y  -  a 

(49) 
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The  transformed  solutions  of  the  fields  should  also  satisfy  (53) 
and  (56),  that  is 


Ez(x*y,kz)  +  Ez  (x»y»kz)  +  Ez  .)  2  2  2  "  0 


(57) 


x  +  y  «  a 


■pri/ 


E  (x.y.k  )  +  E  (x,y,k  )  +  E  (x,y,k  )  - 

z  z  z.  z  z,  z  /  ,  ,  ,  x  2  2 

1  2  x  +  (y  -  d)  «  a 


0  (58) 


Substituting  the  appropriate  transformations  in  a  manner  analogous  to 
that  in  "Formulation  of  the  Problem,"  (57)  and  (58)  may  be  cast  in  a 
matrix  form  similar  to  that  which  occurs  in  the  familiar  method  of 
moments  formulation  [3].  We  have 


p 

* 

- 

V 

2 12 

J1 

V1 

_Z12 

Z22_ 

- 1 

CN 

-3 

_ 1 

II 

1 - 

< 

to 

(59) 


where  the  elements  of  the  matrix  [Z ]  are  called  "generalized  impedances 
and  are  given  by, 


2  2 
k  -  k 


'll 


- - 51  H<2)  (/k2"-  k2  a),  Z„  =  Z 

4  (i)i.  o  z  22  11 


k2  -  k2 


(60) 


X, 


'■  H(2)  -  k* 


d) 


‘12  4  oif.  o 

t lie*  elements  of  the  vector  £  are  called  "generalized  currents"  and  are 
given  by 

VVV’VVV  (61) 

and  the  elements  of  the  vector  $  are  called  "generalized  voltages"  and 
are  given  by 

(62) 


rA 


1 


1 


*  V2  =  “Ez. 
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In  order  Co  determine  Che  modal  solution  of  the  problem,  which  is 
connected  with  the  zeros  of  the  total  impedance  of  the  system,  we  seek 
for  all  possible  source-free  solutions  of  (59).  Hence  t  ■  appropriate 
matrix  equation  is  (59)  with  the  right  side  zero;  that  is 


Z11 

* 

Z12 

V - 

1— 1 

1 _ 

•  ■ 

0 

Z]2 

zn_ 

A 

0 

•  « 

Equation  (63)  will  have  a  non-trivial  solution  provided  that 
the  impedance  matrix  is  zero.  That  is 


(63) 

the  determinant  of 


which  may  be  rewritten  as 


(64) 


Z11  *  Z12 


(65) 


Substituting  (60)  into  (65)  results  in 


k  =  ±k  (66) 

z 

which  represents  the  TEM  mode;  and, 

H(2)  (A2  -  k'2  a)  +  H(2)  (A2  -  k2  d)  -  0  (67> 

o  z  o  z 

which  is  identical  to  (38).  Once  the  proper  values  of  k  for  which  (63) 
is  satisfied  are  found,  the  eigenvectors  J^,  Jj  may  be  obtained  by  solving 
(59)  for  each  k . 
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VII.  Results  and  Conclusions 

The  case  of  plane  wave  excitation  is  shown  in  Figures  7-9. 

The  example  treats  an  incident  plane  wave  for  which  90  ,  the  polar 
angle  is  45°  and  <t>  ,  the  azimuthal  angle  is  45°.  Figures  7-8  describe 
the  magnitude  of  the  currents  induced  on  the  wires  -  Fqs.  (35), (36)  - 
vs.  the  normalized  distance  between  the  wires,  for  a  fixed  ratio 
d/a  -  10.  Note  the  oscillations,  which  reflect  the  changes  that 
occur  in  the  mutual  coupling  between  the  wires  as  we  vary  the  distance. 
Figure  9  shows  the  magnitude  of  the  currents  induced  on  the  wires 
vs.  the  ratio  d/a  for  a  given  distance  between  the  wires  d/A  *•  2.5. 

One  may  observe  that  the  currents  decrease  monotonically  as  the 
ratio  d/a  is  increased  (i.e.,  the  radius  a/A  is  decreased). 

The  case  of  voltage  source  excitation  is  shown  in  Figures  10 
through  20.  Figures  10-13,  show  the  variation  of  the  imaginary  part 
of  the  normalized  z-directed  wave  numbers,  which  are  the  normalized 
attenuation  constants  in  the  z-direction,  as  a  function  of  the  normalized 
distance  between  the  wires,  for  a  given  ratio  d/a.  We  have  depicted 
only  the  two  lowest  even  modes  and  the  two  lowest  odd  modes,  and 
we  have  chosen  d/a  *  10  and  d/a  *>  100  for  Figures  10-11,  and  Figures  12-13 
respectively.  By  looking  at  the  graphs  one  can  readily  see,  as  may 
be  expected,  that  for  fixed  d/a,  the  attenuation  constants  decrease 
as  the  distance  between  the  wires  is  increased.  In  contrast  the 
attenuation  constants  increase  for  fixed  distance  d/A  as  the  ratio 
d/a  is  increased  (i.e.  when  the  radius,  a/A,  is  decreased).  This 
fact  is  also  demonstrated  in  Figure  14  which  shows  the  variation  of 
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le  currents  Induced  on  the  wires  due  to  plane  wave 
:citation,  vs.  the  normalized  distance  between 
ie  wires  (d/a  =  10,  d/A  <  1) • 


i  tne  vires  due  to  plane  wave 
mail zed  distance  between 


The  currents  induced  on  the  wires  doe  to  plane  wave 
excitation  vs.  the  ratio  d/a*  (d/X  *  2.5). 
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Fig.  13  The  normalized  attenuation 
even  and  odd  modes  on  the  \ 
distance  between  the  wires 


ODO  MODE 


Fig.  14  The  normalized  attenuation  constants  of  the 

lowest  even  and  odd  modes  on  the  wires  vs.  the 
ratio  d/a,  (d/X  *  2.0). 
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The  lowest  odd-mode  current  at  z  *  0  normalized 
to  the  TEM  current  vs.  the  normalized  distance' 
between  the  wires  (d/a  *  10,100;  d/X  <  1). 


lowest  odd-mode  current  at  z  *  0  normalized 
he  TEM  current  vs.  the  normalized  distance 
reen  the  wires  (d/a  *  10,100;  d/X  >  1). 
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Fig.  Al.  Illustration  of  the  contours  in  the  k  plane. 
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Also,  hJ  Is  ths  Hsnkel  function  of  the  second  kind  of  first  order,  and 


we  have  used  the  relation 


"£2><x)  *  -»J2> 


Using  the  residue  theorem  the  integral  in  Equation  (Al)  results  in  the 
residue  at  the  simple  pole  k  •  k  ..  thus  we  have, 

.4,  k 


,  .  11  k  _  -  _ 

P‘  "  ^  -  k2p,  V  !(■  h<2)(,42  -  k2  .)  1  J  «'2)(A2  -  k25id)] 


Ik  ,  t 


where 


n  •  <^)1/2 


Substituting  k  (*  //?  -  k^p+)  into  (A5)  one  readily  obtains  Equation  (40). 


A2. 


The  Contribution  of  the  Branch  Point 


In  this  appendix  the  contribution  of  the  branch  point,  namely 
Equation  (41),  is  evaluated. 

The  Integral  around  the  branch  point,  as  shown  in  Figure  Al,  is 
given  by, 


rCDE 


-1 


¥  3“ 


2  2  2 
x  +  y  -  a 


CDE 


(k2  -  k2)  [H^2>(v42  -  a)  ±  H<2)(*<5  -  k2  d)] 

2  0  Z  O  Z 


jk  t 


dk 


By  changing  the  variable  of  integration  to 


k  »  -k  +  ee 

z 


je 


where  6  varies  from  -  ~V  to  we  obtain 


(A6) 


(A7) 


k2  _  y2  m  (k  +  k  )  (k  -  k  )  -  ee36  (2k  -  eej6) 
z  z  z 


(A8) 


In  the  limit  where  e  -*■  0  we  have. 


E*nc  (x,y,kz) 


2  2  2  ^inc 

:c  +  y  »  a 


*-  ,  a  4. 

x  +  y  ■  a 


(A9) 


k2  -  k2  -*•  2kee36 

Z 


jr^? 


Avz  e3® 


H^2)(/k2  -  k2  a)  ■>  1  -  j-  Hn 
n  z  ^ 


H(23(/k2  -  k2  d)  ->•  1 


(A10) 

(All) 

(A12) 

(A13) 


,DEF 


Substituting  Equations  (A9-A13)  in  Equation  (6),  we  obtain  I+  ‘  ■  0 
i.e.,  zero  contribution  for  symmetrical  excitation,  and, 


(A14) 
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,DEF  _  2ui€  -pine 


Iw  -  ~  E*  (x,v,-k)  ?  -  . 

z  x*  +  y*  -  a* 


2kc  #je  (+jj  *n  J> 


je  *j0d6 


n  £n( 


-t—  ljnc  (x,y,-k)  2  2 

f)  *  x2  +  y2  «  a2 


(Ali) 


where 


"  -  <*>'» 


It  is  easily  seen  that  Equation  (A15)  can  be  written  in  the  explict  form. 


l“F  -  <Zo)"1  E*nc  (x,y,-k) 


2,2  2 


n  d 

Where  Zq  *  Cn  is  the  characteristic  impedance  of  the  two-wire  line 
as  given  by  ( 9 ] , 


x  +  y  ■  a 


By  changing  cha variable  of  integration  in  tha  same  way  aa  Equation 
(A18)  whara  u  variaa  from  R  to  t  ,  ona  obtaina 


-  u2(l  +  j£) 


XT 


-u(l  +4)1/2 
u 


4* 


(A23) 


(A24) 


V.  A  ►  k 

where  9  ■  tan  — 

Taking  tha  limit  whara  £  0  and  R  »  equation  (A22)  results  in 


,BC 


lim 


*  0 
-►  oo 


E*nC  (x,y ,-k  +  ju) 


-jkz  -u* 

_ lx2  -k  y2  ,  a2 

u2  (1  +  i~)  (H^2)[-u(1  +  ^-)1/2  ej<<>/2  a]  ±  H^2)[-u(l  +  ^|)1/2  eJ<))/2 


(A25) 
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Appendix  B 

The  Method  for  Finding  the  Roots  of  Equation  (38) 

A  well  known  and  proven  method  for  finding  the  roots  of  equation 
(38),  is  that  of  the  steepest  decent.  Equation  (38)  may  be  modified 
using  a  real  function  h  in  the  following  fashion, 

<ki>  - 1 W  I  - 0  <B1> 


where  h  (k  )  is  given  by  Equation  (A2)  and  k  the  z-directed  axial  wave 
i  z  z 

number  may  be  cast  into  the  vector  form 


k 

z 


(B2) 


t  z.t 

where  ft  and  tj.  are  the  real  and  imaginary  parts  of  k  ,  respectively, 
z  z  z 

The  search  by  this  method  starts  with  initial  value  1;  for  the 

zo 

wave  number  vector.  The  normalized  gradient  of  the  absolute  value  of 
h  (k  ),  obtained  by  differentiating  (Bl) : 

±  7. 


3habs+(kz> 
Sign  (  3ft  * 


'dh  ,  (k  ) 
abs+  z 

Sign  ( - - - ) 


(B3) 


is  measured  and  the  vector  k  is  altered,  in  accordance  with  the  negative 

z 

of  the  value  obtained.  This  procedure  is  repeated  causing  the  value  of 

h  to  be  successively  reduced  and  the  t.  vector  to  approach  the  root 
abs,  z 

of  Equation  (Bl). 

(■'or  convenience,  the  +  and  -  subscripts  will  be  omitted  in  the 
following  discussion.  The  method  can  be  described  by  the  relation 


t 

\  +  1 


t 


[uk]\ 


(BA) 


A7 


where 


8h  .  (k  ) 
Sign  > 


Sign  ( 


*habs<V 


k  -  k 

2  Z, 


Sign  (z) 


1  x  vi 

m 

-1  x  <  i 

» 


is  the  normalized  gradient  at  a  point  in  the  k^  plane  corresponding  to 


k  •  k  and 
z  zk 


3/4  +  i/4  Vj.j  ?£ 


pk  -  1 


3/4  +  1/4  1  2 

k  -  1  k 


is  a  matrix  that  together  with  an  initial  matrix [^J control  stability 
and  rate  of  convergence. 

This  method  is  suitable  for  our  purpose,  since  we  are  interested 

only  in  leaky  modes  having  the  smallest  attenuation  in  the  z  direction. 

Thus,  we  are  allowed  to  sea  ch  for  the  zeros  within  a  rectangular 

region  of  the  second  quadrant  in  the  complex  k  (  *  8  +  id  )  plane, 

z  z  z 

which  is  defined  by 


0  <  <  k 


0  <  a  «  k 

z 


The  absolute  value  of  the  complex  function  Equation  (Bl)  is  reduced  to 


less  than  10  at  the  location  of  a  root. 


48 


The  first  two  programs  calculate  the  induced  currents  in  the 
case  of  a  plane  wave  excitation.  The  first  program  varies  the 
distance  between  the  wires  d  (in  wavelengths)  while  keeping  the 
ratio  d/a  constant.  The  second  program  varies  the  ratio  d/a  while 
keeping  the  distance  between  the  wires  constant. 

The  following  four  programs  find  the  real  and  imaginary 
parts  of  the  wave  numbers  associated  with  higher  order  modes  in 
the  case  of  volt  excitation,  and  calculate  the  contribution  of 
each  mode  to  the  total  current.  The  first  program  varies  the  distance 
between  the  wires,  while  keeping  the  ratio  d/a  constant,  for  the 
even  leaky  modes  and  the  second  carries  it  out  for  the  odd  leaky 
modus.  The  third  program  varies  the  ratio  d/a,  while  keeping  the 
distance  between  the  wires  constant ,  for  the  even  leaky  modes  and 
the  fourth  carries  it  out  for  the  odd  modes. 


I 


l 

t 

i 

t," 


I 

h 

l! 


c  . ... . . . . 

C  THE  PURPOSE  OF  THIS  PROGRAM  IS  TO  CALCULATE  THE 

C  INDUCED  CURRENTS  IN  THE  CASE  OF  A  PLANE  WAVE 

C  EXCITATION. 

C  THE  PROGRAM  VARIES  THE  DISTANCE  BETWEEN  THE  WIRES- 

C  * D " ( I N  WAVELENGTHS)  WHILE  KEEPING  THE  RATIO  "D/A* 

C  CONSTANT 

C  . . . . 

COMPLEX  KAf KD»CQ»CE*CUi 
COMPLEX  HOA » HOD  >  JOA r JOD  ? YOA r YOD 
COMPLEX  CU2 

REAL  K» I JOA» IYOAf I  JOD* I YOD 
C 

C  THE  RATIO  •D/A' 

C 

M-100 

WRITE<5»50)  M 

50  FORMAT <///* IX* '  D/A  IS  '»I4»///> 

PI«3 ♦ 1415927 
K*2.*PI 
C 

C  THE  VERTICAL  AND  AZIMUTHAL  ANGLES  OF  THE 

C  INCIDENT  WAVE 

C 

Q-PI/4 

PY-PI/4 

C 

C  INITIAL  VALUE  FOR  THE  DISTANCE-!) 

C 

D-0. 1 
10  A=D/M 


l 


40 


C 

c 

c 

c 


IF(D.GT.S.OO)  GO  TO  20 

WRITE  (5*40)  D 

FORMAT (/flX* '  D  IS  ' rP6. 4) 

KA=K*A*SIN<a)*< 1 ♦ »0< ) 

KD»KA*M 

CALL  CBESO  <  KA  » JOA, YOA) 

CALL  CBESO (KD  » JOD, YOD) 

R JOA-REAL  <  JOA ) 

IJOA*AIMAG< JOA) 

RYOA-REAL(YOA) 

iyoa=a;i:mag(yoa> 

R  JOD-REAL  ( ,)0D  > 

I JOD- A I MAGI JOD) 

RYOD=REAL ( YOD) 

I YOD-AIMAG  <  YOD ) 

THE  RELATIVE  CONTRIBUTION  OF  THE  EVEN  AND  ODD  PARTS 
AND  THE  TOTAL  CURRENT  INDUCED  ON  THt  WIRES 


■* 

i 


I  i 

i  * 

I  \ 
!  i 


HOA=CMPLX <  RJOA+ 1 YOA  *  I JOA -  R YOA ) 

HOD=CMPLX  <  RJOD+IYOD  r I JOD--RYOD) 

CE=2 ♦ /377 . /K/SIN < 0) #  ( 1+CEXP ( < 0 . .  1. .  >*KD*SIN<PY> )  )/<HOAFMOD> 
C0*2./377./K/SIN<Q)*  <  1--CEXP  ( <  0  < . » 1 .  >*KD*8IN(PY> ) )  /  ( HO  A-  HOD ) 


CU1=CE+C0 
ACU1-CABS(CU1 ) 
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RAT-CABS (CE/CO) 

WRITE<5»30) ACUlrRAT 

30  FORMAT< IX > 'CURRENT  1  IS'  »E15.8>/»1X» 'EVEN/ODD  IS'»E15.8 

CU2-CE-C0 
ACU2-CABS ( CU2 ) 

WR I TE  <  5 » 37 )  ACU2 

37  FORMAT < IX > 'CURRENT  2  IS',E15.8) 

RAT1-ACU1/ACU2 

WRITE(5» 38)  RATI 

38  FORMAT < lXi ' CUR1/CUR2  IS' »E15.8»///) 

C 

C  CHANGING  THE  DISTANCE-D 

C 

D-D+O.l 
GO  TO  10 
20  CONTINUE 

STOP 
END 
C 

C  THE  SUBROUTINE  CBESO  CALCULATES  THE  ZERO  ORDER  BESSEL 

C  FUNCTIONS  FOR  A  GIVEN  INPUT  VARIABLE  Z 

C 

SUBROUTINE  CBESO < Z,BSJO» BSYO > 

COMPLEX  Z » BS JO , BSYO » Y » W » PO » FO * CEXP 

COMPLEX  CCOS  nCSINn CSQRT , CLOG , SPO » SQO , QO * CS » SN » SP 1  *  SO 1 

PI=3. 141593 

A=CABS<Z) 

BSJ0-C1 .  »0. ) 

IF(A.EQ.O.)  GO  TO  2 
IF< A.GT.3. )  GO  TO  1 
Y=Z*Z/9. 

BS J0=1 ♦ +Y* ( -2 . 2499997+Y* ( 1 ♦ 2656208+Y* ( - . 3163866 
1+Y*< .0444479+Y*<~. 0039444+Y#. 00021) ) ) ) ) 

GO  TO  2 

1  U=3./Z 

F0= ♦ 79788456+W*  < - ♦ 00000077+W#  < - ♦ 0055274+W# ( - . 00009512 
1+W*  < . 00137237+W* ( - . 00072805+W* . 00014476 ) > ) ) ) 

P0=  *785398 16+W*  < . 04166397+W# < ,00003954+W*<-. 00262573 
1+U*<  *00054125+W*<  .00029333--U*. 00013558) ) ) )  > 

BS JO=FO*CCOS  <  Z-PO ) /CSQRT ( Z  > 

2  CONTINUE 
BSY0=<-i.E30*0) 

IF<A.EQ.O.)  GO  TO  3 
IF  <  A ♦ GT ♦ 3 ♦ )  GO  TO  4 

BSYO- .  63661977*CL.0G <  *5*Z )  #BS J0+ . 36746691+Y*  ( .  60559366 
l+Y*<-.74350384+Y*<  *2530011?+Y*<--.04261214+Y*< .00427916 
2-Y*.00024846>)))) 

GO  TO  3 

4  CONTINUE 

BSY0=F0*CSIN ( Z-PO ) /CSQRT  <  Z ) 

3  CONTINUE 
RETURN 
END 
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Output  Samples  for  Program  No 


D/A  IS  10 


D  IS  .1000 

CURRENT  1  IS  0.62013944E-03 
EVEN/ODD  IS  0.29726247E+01 
CURRENT  2  IS  0. 89227682E-03 
CUR1/CUR2  IS  0 . 6950079 1E+00 


D  IS  . 2000 

CURRENT  1  IS  0.71889397E--03 
EVEN/ODD  IS  0 . 20512333E+01 
CURRENT  2  IS  0. 12857444E-02 
CUR 1/ CURS  IS  0. 55912665E+00 


1)  IS  .3000 

CURRENT  1  IS  0.81891406E--03 
EVEN/ODD  IS  0.17458486E+01 
CURRENT  2  IS  0 . 162156 14E--02 
CUR1/CUR2  IS  0. 50501 576E+00 


D  IS  .4000 

CURRENT  1  IS  0 ♦  9581 0231 E--03 
EVEN/ODD  IS  0.15963233E+01 
CURRENT  2  IS  0 . 19070078E-02 
CUR 1 /CURS  IS  0. 50241 131E+00 

I 


D  IS  .5000 

CURRENT  1  IS  0.11S98827E-02 
EVEN/ODD  IS  0. 14861194E+01 
CURRENT  2  IS  0. 21 185441E-02 
CUri/CUr2  IS  0.54749046ET00 


D  IS  .6000 
CURRENT  1  is 
EVEN/ ODD  IS 
CURRENT  2  IS 
CUR 1/ CURS  IS 


0.14269433E-02 
0. 13361567ET01 
0.217A4565E-02 
0 . 655&2683E+00 
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nnn  noon  non  ooooooo 


program  No.  2 


THE  PURPOSE  OF  THIS  PROGRAM  IS  TO  CALCULATE  THE 
INDUCED  CURRENTS  IN  THE  CASE  OF  A  PLANE  WAVE 
EXCITATION. 

THE  PROGRAM  VARIES  THE  RATIO--D/A* »WHILE  KEEPING 
THE  DISTANCE  BETWEEN  THE  WIRES-*D*  C0N8TANT 

COMPLEX *KA»KD> CO »CE»CU1 ******* 

COMPLEX  HOA  » HOD  t JOA  *  JOD » YOA » YOD 
COMPLEX  CU2 

REAL  K» IJOA» IYOA» I  JOD* I YOD 
INITIAL  VALUE  FOR  THE  RATIO  'D/A* 


M=10 

PI=3. 1415927 
K--2.*PI 

THE  VERTICAL  AND  AZIMUTHAL  ANGLES  OF  THE 
INCIDENT  WAVE 


G=PI/4 
PY=PI/4 

THE  DISTANCE  BETWEEN  THE  WIRES 
D=2 .5 

WRITE  <5»^0)  D 
40  FORMAT </* IX* 'D  IS  ' *F4.4) 

10  A=D/M 

IF(M.GT.IOOO)  GO  TO  20 
WRITE<5*50)  M 

50  FORMAT (///» IX* '  D/A  IS  '*14) 

KA=K*A*SIN<G)*<1.*0.) 

KD=KA#M 

CALL  CBESO(KA  »JOA»YOA) 

CALL  CBESO(KD  * JOD* YOD) 

R JOA^REAL  <  JOA ) 

IJOA=AIMAG< JOA) 

RYOA=REAL< YOA) 

IYOA=AIMAG< YOA) 

RJ00=REAL< JOD) 

I JOD=AIMAG< JOD) 

RYOD=REAL< YOD) 

IYOD=AIMAG< YOD) 

HOA=CMPLX<RJOA+IYOA* IJOA-RYOA) 

HOD=CMPLX<RJOD+IYOD* IJOD-RYOD) 

C 

C  THE  RELATIVE  CONTRIBUTION  OF  THE  EVEN  AND  ODD  PARTS 

C  AND  THE  TOTAL  CURRENT  INDUCED  ON  THE  WIRES 

C 

CE-=2./377./K/SIN<G)#  ( 1+CEXP<  (0. * 1 .  >*KD*SIN<PY) )  >/<HOA+HOD) 
C0=2 ♦ /377 . /K/SIN ( G ) #  <1-CEXP< (0. * 1 . )*KD#SIN<PY> ) )/(HOA-HOD) 
ClJl-CE+CO 
ACU1-CABS<CU1 ) 

RAT=CABS(CE/CO) 
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M=M+1 
GO  TO  10 
20  CONTINUE 

STOP 
END 
C 

C  THE  SUBROUTINE  CBESO  CACULATES  THE  ZERO  ORDER  BESSEL 

C  FUNCTIONS  FOR  A  GIVEN  INPUT  VARIABLE  Z. 

C 

SUBROUTINE  CBESO<Z»BSJO»BSYO) 

COMPLEX  Z » BS JO , BSYO » Y » U . PO * FO » CEXP 

COMPLEX  CCOS , CSIN » CSGRT *  CLOG  * SPO  t SOO » QO , CS » SN » SP1 » SQ1 

PI=3. 1.41593 

A=CABS<Z) 

BSJO=< 1 « »0. ) 

IF< A.EQ.O. >  GO  TO  2 
IF<A.GT.3.)  GO  TO  1 
Y=Z*Z/9. 

BSJ0=1 . +Y* ( -2 . 2499997±Y*  < 1 . 2656208+Y*  < - .3163866 
1 +Y*  < . 0444479+Y*  < - . 0039444+Y* . 00021 ) ) > ) ) 

GO  TO  2 

1  Wa3./Z 

F0« . 79788456+W#  < - . 00000077+W*  < - . 0055274+U*  < - . 00009S1 2 
1+W*< .00137237+W*<-.00072805+W*. 00014476) >  > ) ) 

PQ= . 78539816+W* < .04166397+W#( .00003954+W*<- .00262573 
1+W*< .00054125+W*< .00029333-W*. 00013558))))) 

BS J0=F0*CC0S ( Z-PO ) /CSQRT <  Z ) 

2  CONTINUE 
BSY0-<-i .E30»0> 

IF (A.EQ.O. )  GO  TO  3 
IF(A.GT.3.)  GO  TO  4 

BSYO85 . 63661977#CL0G  < .  5*Z )  *BS J0+ .  36746691+Y*  <  .  60559366 
l+Y*<-.74350384+Y*< .25300117+Y#<-.04261214+Y*< .00427916 
2- Y* .00024846) ) ) ) ) 

GO  TO  3 

4  CONTINUE 

BSY0=F0#CS IN  < Z-PO ) /CSQRT  <  Z ) 

3  CONTINUE 
RETURN 
END 


3 


Vi 

i 

■'i 
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Output  Saaple  for  Progrta  No.  2 


D  IS  2*5000 


D/A  IS  10 

CURRENT  1  IS  0.40766686E-02 
EVEN/ODD  IS  0*16611102E'tQl 
CURRENT  2  IS  0.28035670E-  02 
CUR1/CUR2  IS  G«  14541007E t 01 


D/A  IS  11 

CURRENT  1  IS  0.38489985E  02 
EUEN/ODD  IS  0.1547599SET01 
CURRENT  2  IS  0.2S430033E-02 
CUR1/CUR2  IS  0»1S130081E+01 


D/A  IS  12 

CURRENT  1  IS  0 . 36601 04 IE'  02 
EVEN /ODD  IS  0 . 146256901. 101 
CURRENT  2  IS  0,235610731-02 
CUR1/CUR2  IS  0 . 155340112+01 


D/A  IS  13 

CURRENT  1  IS  0.35C0931TE -02 
EVEN/oDD  IS  0.13933734E+01 
CURRENT  2  IS  0 .22I023I7E-C2 
CUR 1 /CUR2  I 5  0 . 15796775ET01 


D/A  IS  14 

CURRENT  1  IS  0 ♦ 336464 15E -02 
EVEN/ODD  IS  0 »  J337432GE101 
CURREN \  2  IS  0.210S6C1SE-02 
CCR1/CUR2  IS  0 • 1 59561 37ETC1 
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Progrui  Mo.  3 


C  . . . . . . . . . . . 

C  THE  PURPOSES  OF  THIS  PROGRAM  ARE! 

C  1)  TO  FIND  THE  REAL  AND  IMAGINARY  PARTS  OF  THE  REDIRECTED 

C  AND  Z-DIRECTED  WAVE  NUMBERS .ASSOCIATED  WITH  THE  HIGHER 

C  ORDER  EVEN  LEAKY  MODES  IN  THE  CASE  OF  VOLT  EXCITATION 

C  2>  TO  CALCULATE  THE  CONTRIBUTION  OF  EACH  MODE  TO  THE  TOTAL 

C  CURRENT. 

C  THE  PROGRAM  VARIES  THE  DISTANCE  BETWEEN  THE  WIRES- 

C  *D*-< IN  WAVELENGTHS)  WHILE  KEEPING  THE  RATIO  ■D/A* 

C  CONSTANT 

C  . . . . . . . . . . . 

COMMON  KR.AA.M.BR.AR 
COMPLEX  KR.KZP.KRP.RES 

REAL  I JOA  » I JOB. IJ1A. I J1D.IYOA.IYOD.IV1A. IY1D 
COMPLEX  KRA » KRD » JOA » YOA  » JOD  » YOD  »  HOA » HOD » J1 A 

COMPLEX  Y1A.J1D.Y1D.GR1.GR2.H1A.H1D.F.GRA1.GRA2 
COMPLEX  KRPA.KRPD 
REAL  KZ<2) »K»KH 
REAL  GRAD<2) 

REAL  OGRAD <  2 ) . CON ( 2 ) 

C 

C  DEPARTURE  POINT  FOR  THE  SEARCHING  IN  THE  KZ  PLANE 

C 

1*1 
BR*1 » 

AR*-2 . 

C 

C  START  SEARCHING  FOR  THE  ROOTS  IN  THE  KZ  PLANE 

C 

9  IF  (BR.GT.lOO)  GO  TO  8 

WRITE  <5»18)  I 

18  FORMAT <.///.  1 X  »'  THE  INDEX  OF  THE  MODE  IS  '.12.//) 

OGRAD < 1 )*0. 

0GRAD<2)=0. 

CON( 1 )=0,50 
C0N<2)*0.50 
N*2 

PI-3.1415927 

K»2*PI 

KH«K##2+AR**2-BR**2 

BZ*SQRT<0.5*<KH+SQRT<  <KH**2+<2*AR*BR>**2> ) ) ) 

AZ— AR#BR/BZ 
KR-CMPLX  <  BR » -AR ) 

KZ< 1 )*BZ 
KZ<2)»AZ 
C 

C  THE  RATIO  *D/A* 

C 

M*100 

AA=1./M 

WRITE  <5.19)  M 

19  FORMAT  <1X»'  D/A  =  '»I6  ) 

3  CALL  FUNCT <N»NZ»VAL.GRAD> 

IF(VAL.LT .0. IE-5)  GO  TO  1 
IF  (GRAD( 1 > )  22.27.27 
22  IF  < OGRAD < 1 ) )  23.24.24 
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CON ( 1 ) »CON ( 1 > /2 • 

KZ ( 1 >  «KZ  < 1 >  +CON ( 1 ) 

OCR AD ( 1 >  «GR  AD < 1 > 

00  TO  31 

IF  ( OQRAD < 1 ) >  28*29*29 
C0N(l)-C0N<l>/2 
KZ  < 1 >  »KZ  < 1 > -CON  < 1 > 
OQRAD  < 1 >  *GRAD  < 1 ) 
CONTINUE 

IF  <0RAD<2>>  32*37*37 
IF  (0QRAD<2) )  33*34*34 
C0N(2>»C0N<2)/2. 

KZ  <  2 ) *KZ  <  2 ) +CON  <  2  > 
OQRAD  <  2 ) *GRAD ( 2  > 

GO  TO  3  ' 

IF  <0GRAD<2) )  38*39*39 
C0N(2>=C0N<2>/2 
KZ ( 2 ) “KZ  <  2 )  ~CON  ( 2  ) 
OQRAD  <  2 ) “GRAD  <  2  > 

60  TO  3 
CONTINUE 


INITIAL  VALUE  FOR  THE  DISTANCE-D 


D=Q.  1 
A=D/M 

IF<D»GT.5«00)  GO  TO  60 

BP=+BR/B 

AP=+AR/D 

KRP=CMPLX  <  BP  » -AP ) 

NRPA=NRP*A 

KRPD=KRP*D 

K=?*PI 

KH«K**2+AP**2-BP**2 

BZ=SQRT <  0 . 5*  < KH+SQRT < ( KH**2+  <  2*AP*BP ) #*2  >  > ) ) 

AZ=-AP*BP/BZ 

KZP=CMPLX<-BZ*AZ) 

WRITE  (5*80)  D 

FORMAT < //* IX* 'D  IS  *  * »F6 *4) 

URITE<5*90)  KZP 

F0RMAT</*1X» 'REAL  KZ  IS  t'*F9.4*'  I MAG  KZ  IS  i'*F9»4) 

CALL  CBESKKRPA* J1A*Y1A> 

CALL  CBES1 (NRPD* J1D»Y1D> 

RJ1A=REAL< J1A) 

I J1A»AIMAG< J1A) 

RY1D=REAL< Y1D) 

TY1D«AIMAG< Y1D) 

RJ1D«REAL( J1D> 

IJ1D=AIMAG( J1D) 

RYlDaREAL(YlD) 

IY1D=AIMAG(Y1D) 

H1A=CMPLX<RJ1A+IY1A* I J1A-RY1A) 

HU«=CM*LX<RJ1D+IY1D*  +  IJ1D-RY1D> 

THE  CONTRIBUTION  OF  THE  LEAKY  MODE 
TO  THE  TOTAL  CURRENT  NORMALIZED  TO  THE 


u  uu  >oooou  as  uuoo 


C  CONTRIBUTION  OF  THE  TEH  MODE 

C 

RES»4./377. 

RES*  RES*<0.1)*K/KRP/KZP/<A#H1A+D*H1D> 

ARS=CABS<RES> 

T*M 

RAT  *ARS#377 ♦ #ALQG <  T  > /P I 
WRITE <5 .40)  ARS 

40  FQRMATI/r 1X» 'THE  CURRENT  IS  t'»El5.8) 

WRITEI5.41)  RAT 

41  FORMAT </. 1X» 'CUR/TEM  ISI'.EIS.B) 

CHANQINO  THE  DISTANCE-D 

D*D+0» 1 
GO  TO  70 
0  CONTINUE 

CHANGING  THE  DEPARTURE  POINT  IN  ORDER  TO 
START  THE  SEARCHING  FOR  THE  CONSECUTIVE  ROOT 

1*1  +  1 
BR=BR+5 
GO  TO  9 
CONTINUE 
STOP 
END 

THE  SUBROUTINE  CBESO  CALCULATES  THE  ZERO  ORDER  BESSEL 
FUNCTIONS  FOR  A  GIVEN  INPUT  VARIABLE  Z 

SUBROUTINE  CBES0<Z,BSJ0,BSY0> 

COMPLEX  Z»BSJO*BSYO» Y»W»PO»FO»CEXP 

COMPLEX  CCOS » CSI N » CSQRT » CLOG , SPO , SQO , 00 » CS r SN » SP 1 » SQ1 

PI=3. 141593 

A»CABS(Z> 

BSJ0*(1*  »0. ) 

IFIA.EQ.O.)  GO  TO  2 
I F  <  A  ♦  GT  ♦  3 . )  GO  TO  1 
Y*Z*Z/9. 

BSJ0*1 . +Y*  < -2 . 2499997+Y*  < 1 . 2656208+Y* ( - . 3163866 
1+Y#< .0444479+Y*< - .0039444+Y# . 00021 ) ) > ) ) 

GO  TO  2 

1  (4*3 . /Z 

FO* .  797884!:  6+W#  <  - .  00000077+W*  ( - .  0055274+W#  ( - .  00009512 
1+U*< . 00137237+W*<-.00072805+W#. 00014476) ))) ) 
P0*,7853?816+U#<  .  04166397+W*< . 00003954+W*(- . 00262573. 
1+W*< . 00054125+W#  < .00029333-W*. 00013558) ) ) ) )  k 

BS JO*FO*CCOS  <  Z  -PO )  /CSQRT  ( Z ) 

2  CONTINUE 
B8Y0*<-1 *E30» 0) 

IF  <  A ♦ EQ . 0 ♦ )  GO  TO  3 
IF  (A.  GT.  3.)  (30  TO  4 

BSY0*.63661977*CL0G< .5*Z)*BSJ0+.36746691+Y*< .60559366 
1+Y*<~.74350384+Y*< .25300117+Y*<-.04261214+Y*< .00427916 
2--Y* ♦ 00024846 ) ) ) ) ) 


i 


\ 


•  I 
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00  TO  3 
CONTINUE 

BSYO®FO*CSIN  ( Z  -PO )  /('.SORT  <  Z ) 

CONI INUE 

RETURN 

end 


THE  SUBROUTINE  CBES1  CALCULATES  THE  FIRST  ORDER  BESSEL 
FUNCTIONS  FOR  A  GIVEN  INPUT  VARIABLE  Z 

SUBROUTINE  CBES1 < Z , BSJ1 * BSY1 > 

COMPLEX  Z  t  BSJ.t  »BSY1  *Y*W»P1»F1»CEXP 

COMPLEX  CCDS * CSORT, CLOG »CSIN»SP1» SOI >Gl »F»CS>SNiSP2»SQ2 

PI *3, 141593 

A®CAB8<Z) 

BO J 1*0. 

IF (A.EG.O. )  GO  TO  2 
IF  ( A » GT ♦ 3 ♦ )  GO  TO  1 
Y®Z#Z/9, 

B&J1«Z*<  5+Y*<- ,56249985+Y*< ,21093573+Y*<-,0395  *289 
1  +  Y*< .00443319+Y*<-.00031761+Y*. 00001109)) ) ) ) ) 

GO  TO  2 
W*3,/Z 

FI®*  79788456+W*  < ,00000156+W#< ,01659667+W*< ,00017105 
1+W*<-. 00249511  +W*< ,00113653-W*. 00020033) ) ) ) ) 

PI “2  *  35619449+W*  <  ~ . 12499612+W* ( - , 00005650+W* ( . 00637879 
1+W#  <  -•  •  00074348+W#  < -  •  00079824+ W# ,  00029166 )  > ) ) ) 
BSJ1®F1*CC0S<Z-P1)/CQQRT<Z> 

CONTINUE 

BS Y 1 = ( - 1 • E30 » 0 ) 

IF(A,EQ*0* )  GO  TO  3 
IF(A.GT *3, )  GO  TO  4 

BSYla(-,63661977+Y*< ,2212091+Y*<2. 1682709+Y*<-1 ,3164827 
1+Y*( , 3123951 +Y# ( - , 0400976+Y# . 0027873 ) > ) ) ) >/Z 
2+,63661977*CL0G< ,5*Z)*BSJ1 
GO  TO  3 

BSY1=F1*CSIN<Z~P1)/CSQRT<Z> 

CONTINUE 

RETURN 

END 

THE  SUBROUTINE  FUNCT  CALCULATES  FOR  EACH  VECTOR  NZ 
THE  VALUE  OF  THE  FUNCTIONAL  WHICH  WE  TRY 
TO  MINIMIZE  AND  ITS  GRADIENT 

SUBROUTINE  FUNCT <N*KZ»VAL»GRAD) 

COMMON  KR»AA»M»BR» AR 
DIMENSION  GRAD<  2 ) 

REAL  NZ(2) »K»NH 

REAL  I JOA  , I JOD> I J1A> I J1D» IYOA» IY0D»IY1A»IY1D 
COMPLEX  NR i KRA  * NRD >  JOA » YOA  t JOD » YOD » HOA » HOD » J 1 A 
COMPLEX  Y1A»  J1D* Y1D»  GR1 »GR2»H1A»H1D»F » GRA1 » GRA2 
BZ=KZ< 1 ) 

AZ®KZ(2) 

PI®3. 1415927 
K=2*PI 
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KH=K**2+AZ**2-BZ**2 

BR«SQRT<0.5*<KH+SQRT< <KH#*2+<2*AZ*B2>**2> ) ) > 

AR=-AZ*BZ/BR 

KR=CMPLX  <  BR  *  -AR ) 

KRA=KR*AA 

KRD<=KRA*M 

CALL  CBESO<KRA»  J0A» YOA> 

CALL.  CBESO<KRD>  JODnYOD) 

R JOA-REAL ( JOA ) 

I JOA-AIMAG<  JOA) 

RYOA=REAL<YOA> 

IYOA=A.TMAG<YOA) 

R JOD=REAL ( JOD ) 

I J0D“AIMAG< JOB) 

RYOD=REAL ( YOD ) 

I YOD=AIMAG<  YOD ) 

HOA-CMPLX  <  RJOA+ 1 YOA ,  I  JOA --RYOA ) 

HOD-CMPLX <  RJOD+ 1 YOD . I JOD-RYOD > 

F=HOA+HOD 

VAL..=CABS<F> 

VAL=VAL**2 

CALL  CBESt  <KRAi J1A » Y1A ) 

CALL  CBES1 <KRB» J1D»  Y1D ) 

RJ1A«REAL< J1A) 

I J 1  A5- A I  HAG  ( J 1 A ) 

RY1.D=REAL<  Y1D) 

3'  Y 1  D:"  A I  MAG  ( Y 1 B ) 

RJIDbREAL(JiD) 

IJiD=AIMAG< J1D) 

RY1D*REAL< Y1 D) 

IY1D«AIHAG(Y1D) 

HI  A--CMPLX  <  R J1 A+ 1 Y1 A  »  I J1 A-RY1 A  ) 

H 1 B-CMPLX ( R J 1 D+ 1 Y 1 D » + 1 J 1 D-RY 1 D  > 

GRl--=-Hl  A*  ( “CMPLX  <  BZ  >  ~AZ )  *AA/KR ) 

GRIGORI-HID*  <  -CMPL.X  <  BZ  r  ~AZ )  *AA*M/KR ) 
GR2*GRi*(0*0»~l*0) 

GRA1  =:GR  1  ♦CONFIG  ( F  >  +CQNJG  ( GR 1 )  *F 
GRA2-GR2#C0N JG <  F ) f CON JG ( GR2 ) *F 
GRAD  < 1 ) ~RE AL ( GR A 1 ) 

GRAD  <  2 ) “REAL ( GRA2 ) 

retur: 

END 


Output  Saaple  for  Prograa  No.  3 

D/A  -  10 

D  IS  J  .1000 

REAL  KZ  IS  !  -11.9813  1HAG  KZ  IS  :  32.2049 

THE  CURRENT  IS  J  0 . 409B8265E- 03 
CUR/TEH  ISJ  0. 11325743E+00 

D  IS  :  .2000 

REAL  KZ  IS  J  -6.3111  IMAG  KZ  IS  J  15.2847 
THE  CURRENT  IS  :  0.85170465E-03 
CUR/TEH  ISJ  0.23534023E+00 

D  IS  J  .3000 

REAL  KZ  IS  J  -4.6324  IMAG  KZ  IS  l  9.2551 

THE  CURRENT  IS  J  0 . 13608340E-02 
CUR/TEM  ISJ  0. 37602 118E+00 

D  IS  :  .4000 

REAL  KZ  IS  J  -4.0584  IMAG  KZ  IS  J  5.9422 

THE  CURRENT  IS  J  0 . 1957251 6E-02 
CUR/TEM  ISt  0. 540821 34E+00 

D  IS  J  .5000 

REAL  KZ  IS  J  -4.0681  IMAG  KZ  IS  J  3.7940 

THE  CURRENT  IS  J  0 .25318940E-02 
CUR/TEM  ISJ  0.69960463E+00 

D  IS  J  .6000 

REAL.  KZ  IS  J  -4.4116  IMAG  KZ  IS  J  2.4296 

THE  CURRENT  IS  J  0.27964980E-02 
CUR/TEM  ISJ  0.77271914E+00 


Program  No.  4 


C  . . . . . 

C  THE  PURPOSES  OF  THIS  PROGRAM  ARE* 

C  1 >  TO  FIND  THE  REM.  AND  IMAGINARY  PARTS  OF  THE  REDIRECTED 

C  AND  Z-DIRECTED  NAVE  NUMBERS* ASSOCIATED  NITH  THE  HIGHER 

C  ORDER  ODD  LEAKY  MODES  IN  THE  CASE  OF  VOLT  EXCITATION 

C  2)  TO  CALCULATE  THE  CONTRIBUTION  OF  EACH  MODE  TO  THE  TOTAL 

C  CURRENT. 

C  THE  PROGRAM  VARIES  THE  DISTANCE  BETWEEN  THE  WIRES- 

C  *D'-(IN  WAVELENGTHS)  WHILE  KEEPING  THE  RATIO  *D/A" 

C  CONSTANT 

C  . . . 

COMMON  KR»AA»MrBR*AR 
COMPLEX  KR*KZP*KRP*RES 

REAL  IJOA  * IJOD*I J1A* I J1D* IYOA* IYODr IY1A*IY1D 
COMPLEX  KRA*KRD* JOA*YOA*JOft* YOD*HOA*HOD»JlA 
COMPLEX  Y1A*J1D*Y1D*6R1*GR2*H1A*H1D*F*GRA1*GRA2 
COMPLEX  KRPAt KRPD 
REAL  KZ<2) *K»KH 
REAL  GRAD ( 2 ) 

REAL  OGRAD ( 2 ) *  CON  <  2 ) 

C 

C  DEPARTURE  POINT  FOR  THE  SEARCHING  IN  THE  KZ  PLANE 

C 

1=1 

BR-1. 

AR=-2. 

C 

C  START  SEARCHING  FOR  THE  ROOTS  IN  THE  KZ  PLANE 

C 

9  IF  <BR.GT. 100)  GO  TO  8 

WRITE  <5*18)  I 

18  FORMAT <///* IX »' THE  INDEX  OF  THE  MODE  IS  '*I2*//> 

OGRAD< 1 >=0. 

0GRAD<2)=0. 

C0N< 1 >=0.50 
CON  <  2)=0.50 
N=2 

PI=3. 1415927 
K=2*PI 

KH=K**2+AR*#2-BR**2 

BZ=SQRT  <  0 . 5*  <  KH+SQRT <  <  KH**2+  <  2* AR*BR ) **2 ) ) ) ) 

AZ— ARftBR/BZ 
KR=CMPLX  <  BR  * -AR ) 

KZ< 1 )=BZ 
KZ<2)=AZ 
C 

C  THE  RATIO  *D/A* 

C 

M=100 

AA-l./M 

WRITE  <5*19)  M 

19  FORMAT  <1X*'  D/A  =  '*I6  ) 

3  CALL  FUNCT  <  N » KZ  *  V AL  r  GRAD ) 

IF < VAL *LT .0. IE-5 )  GO  TO  1 
IF  <GRAD< 1 > )  22*27*27 
22  IF  < OGRAD < 1 ) )  23*24*24 
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24  CON< 1 )=CON< 1  )/2» 

23  KZ ( 1 > “KZ ( 1 ) +CQN ( 1 > 

OGRAD < 1 ) “GRAD ( 1 ) 

GO  TO  31 

27  IF  <0GRAD<1>>  20»29»29 

28  CON  < 1 ) “CON  <  1 )  /2 

29  KZ  <  1 )  =KZ  < 1 ) -CON  < 1 ) 

OGRAD ( 1 ) *GRAD ( 1 ) 

31  CONTINUE 

IF  (GRAD(2>)  32»37»37 

32  IF  <0GRAD<2>>  33>34»34 

34  CON ( 2 ) “CON ( 2 ) /2 ♦ 

33  KZ  < 2 ) =KZ ( 2  >  +CON ( 2 ) 

OGRAD  <  2 ) “GRAD ( 2 ) 

GO  TO  3 

37  IF  (0GRA0(2) )  38»39*39  ' 

38  CON  <  2 ) =CON  <  2  > /2 

39  KZ ( 2 ) “KZ  <  2 ) -CON ( 2  > 

OGRAD  <  2 ) “GRAD  <  2 ' 

60  TO  3 

1  CONTINUE 

C 

C  INITIAL  VALUE  FOR  THE  DISTANCE-D 

C 

D=0 « 1 
70  A=D/M 

IF (D.GT *5.00)  GO  TO  60 

BP=+BR/D 

AP=+AR/D 

KRP=CMPLX  <  BP » ~ AP ) 

KRPA“NRP*A 

KRPD=KRP#D 

K=2*PI 

KH==K**2+AP**2-BP**2 

BZ=SQRT <  0 . 5*  < KH+SGRT ( ( KH**2+  <  2*AP*BP  >  **2  >  >  >  > 

AZ=-AP*BP/BZ 

KZP-CMPLX  < -BZ  r AZ ) 

WRITE  <5»80)  D 

80  FORMAT < //» 1X» 'D  IS  :'»F6.4> 

WRITE<5»90)  KZP 

90  FORMAT ( /> 1X» ' REAL  KZ  IS  {',F9.4 r'  IMAG  KZ  IS  J'rF9.4> 

CALL  CBES1 <KRPA» J1A»Y1A) 

CALL  CBES1 <KRPD» J1D» Y1D) 

RJ1 AUREAL < J1A) 

I J1A=AIMAG  <  J1A ) 

RY1D“REAL<Y1D> 

IY1D“AIMAG< Y1D) 

RJ1D=REAL( J1D) 

IJ1D=AIMAG< J1D) 

RY1D=REAL< Y1D) 

IY1D-AIMAG(Y1D) 

HlA“CMPLX(RJlA+IYlAr IJ1A-RY1A) 

H1D=CMPLX < RJ1D+IY1D»  HJ1D-RY1D) 

C 

C  THE  CONTRIBUTION  OR  THE  LEAKY  MODE 

C  TO  THE  TOTAL  CURRENT  NORMALIZED  TO  THE 
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C 


CONTRIBUTION  OF  THE  TEH  MODE 

RES-4./377. 

RE8«  RES#(0»1  >*K/KRP/KZP/<  A*H1A+D*H1D) 
ARS"CABS  <  RES ) 

T=M 

RAT»ARS*377 . *ALOG  <  T ) /PI 
WRITE <5 >40)  ARS 

FORMAT < / » IX  7 'THE  CURRENT  IS  '»E15.8> 
WRITE(5r41 )  RAT 

FORMAT {/ r 1X» 'CUR/TEM  ISJ ' »E15.8> 

CHANGING  THE  DISTANCE-D 

0=0+0 * 1 
GO  TO  70 
CONTINUE 


CHANGING  THE  DEPARTURE  POINT  IN  ORDER  TO 
START  THE  SEARCHING  FOR  THE  CONSECUTIVE  ROOT 

1  =  1  +  1 
BR=BR+5 
GO  TO  9 
CONTINUE 
STOP 
END 


SlEn?V522UIiNE  CBES°  CALCULATES  THE  ZERO  ORDER  BESSEL 
FUNCTIONS  FOR  A  GIVEN  INPUT  VARIABLE  Z 


SUBROUTINE  CBESO<Z*BSJO»BSYO> 

COMPLEX  Z»BSJO»BSYO»Y»W»PO»FO»CEXP 

COMPLEX  CCOS » CSIN  7  CSQRT  *  CLOG  t SPO  r SQO » QO » CS  t SN  t SP1 r SQ1 

PI=3. 141593 

A=CABS(Z) 

B8J0M1,»0.> 

IF<A.EG.O.)  GO  TO  2 
IF<A.GT*3. >  GO  TO  1 
Y=Z*Z/9, 

BS J0=1 ♦ +Y* ( -2 . 2499997+Y*  < 1 , 2656208+Y* ( - . 3163866 
I+Y*( « 0444479+Y#.<-  .0039444+ Y#  .00021 )  )  )  > ) 

GO  TO  2 

1  W=3./Z 

F0= . 79788456+W#  < -  * 00000077+ W4 ( - . 005527 4 +W* ( - , 00009512 
1+U*< . 00137237+W*<-.00072805+W*. 00014476) > ) ) ) 

P0= ♦ 785398 16+W*< .04166397+U*< . 00003954+ W*(~, 00262573 
1+W*< . 00054125+W#< « 00029333-U# . 00013558 ) ) ) ) > 

BS JO=FO*CCOS  <  Z-PO ) /CSQRT  <  Z ) 

2  CONTINUE 
BSY0=<-1.E30»0) 

IF  < A .EQ. 0 . )  GO  TO  3 
IF  < A .GT » 3 ♦ )  GO  TO  4 

BSYO- , 636619774CL0G ( ,5*Z)*BSJ0+.36746691+Y*< .60559366 
1  +  Y#< -  .74350384+Y#< ♦ 25300 1 1 7+Y*(- .04261 21 4+Y*( .00427916 
2-Y*. 00024846 ) ) ) ) ) 
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GO  TO  3 

4  CONTINUE 

BSYO->FO*CSIN<  Z-PO  > /C8QRT <  Z  > 

3  CONTINUE 

RETURN 
END 

THE  SUBROUTINE  CBESl  CALCULATES  THE  FIRST  ORDER  BESSEL 
FUNCTIONS  FOR  A  GIVEN  INPUT  VARIABLE  Z 

SUBROUTINE  CBESl (ZiBSJlxBSYI) 

COMPLEX  Z » BS J1 >  BSY1 » Y » W> PI » F 1 1 CEXP 

COMPLEX  CCOS  f CSQRT  t CLOG » CS I N » SP 1 » SQ1 » Q1 » F » CS  »  SN » SP2 » 8Q2 
Pi»3. 141593 
A»CABS(Z) 

BSJ1=»0. 

IF( A*EQ.O  * )  GO  TO  2 
IF<A.GT.3*>  GO  TO  1 
Y=Z*Z/9. 

BSJ1=Z*( >  5+Y#  < - « 56249985+Y#  < ♦ 21093573+Y4 ( - *  03954289 
1+Y*< . 00443319+Y*<~.00031761+Y*. 00001109) ) ) ) >  > 

GO  TO  2 
W=3./Z 

FI  = « 79788456+W*  < *00000156+W*< *01659667+U*<  *00017105 
1+W*  < -  *  0024951 1+W#  < *001 13653-W#  *00020033 > ) ) ) > 

Pl=2. 35619449+W*  < - . 12499612+W* ( - . 00005650+W* ( . 00637879 
1+W* <-.00074348+W*<-.00079B24+W*. 00029166 ) > ) )  ) 

BS J1=F1*CC0S ( Z-Pl ) /CSQRT  <  Z  > 

CONTINUE 
BSY1® <-l .£30*0) 

IF<A.EQ.O.)  GO  TO  3 
IF<A*GT*3* )  GO  TO  A 

BSY1*  < - . 63661977+Y* ( , 2212091+Y*  <  2 . 1682709+Y*  < -1  *  3164827 
1+Y*< . 3123951+Y*<-.0400976+Y*. 0027873 >  > ) ) >  >/Z 
2+*63661977#CL0G< *5#Z>#BSJ1 
GO  TO  3 

BSY1“F1#CSIN<Z-P1 )/CSQRT  <Z) 

CONTINUE 
RETURN 
END 

THE  SUBROUTINE  FUNCT  CALCULATES  FOR  EACH  VECTOR  KZ 
THE  VALUE  OF  THE  FUNCTIONAL  WHICH  WE  TRY 
TO  MINIMIZE  AND  ITS  GRADIENT 

SUBROUTINE  FUNCT < N ? KZ » VAL » GRAD ) 

COMMON  KR*AA»M*BR»AR 
DIMENSION  GRAD(2) 

REAL  KZ<2) *K»KH 

REAL  IJOA  > IJOD» I JlAx IJIDx IYOAx IYODr IY1A» IY1D 
COMPLEX  KR » KR A f KRD » JOAxYOAx JOD x YOD x HOA x HOD x  J1 A 
COMPLEX  YlAxJIDxYlDrGRI xGR2xH1AxH1DxFxGRA1 xGRA2 
BZ*KZ ( 1 > 

AZ»KZ(2) 

PI=3. 1415927 
K“2*PI 


66 


KH«K**2+AZ**2-BZ**2 

BR-SQRT  <  0 . 5* ( KH+BQRT ( ( KH**2+ ( 2*AZ*BZ ) **2  > ) ) ) 

AR—AZBBZ/BR 

KR»CMPLX  <  BR » -AR ) 

KRA«KR*AA 

KRD-KRABM 

CALL  CBE80 ( KRA » JOA  t YOA  > 

CALL  CBESO  <  KRD  t JOB  »  YOD  ) 

RJOA-REAL(JOA) 

I  JO  A«  A I  MAO  <  JO  A ) 

R YOA=REAL  <  YOA ) 

IYOA»AIMAG<YOA) 

RJOD«REAL( JOD) 

IJOD»AI..AG<  JOD) 

RYOD=*RE^'  i  'YOD) 

IYOD=AIP  )< YOD) 

HOA=CMPLX<  R JOA+IYOA » I JOA-RYOA ) 
HOD=CHPLX(RJOD+IYOD» I JOD-RYOD) 

F=HOA-HOD 

VAL=CABS<F) 

VAL=VAL**2 

CALL  CBES1 (KRA»  J1A» Y1A) 

CALL  CBES1 (KRDf J1D> Y1D) 

RJ1A=REAL( J1A) 

I J1A=AIMAG< J1A) 

RY1D=REAL<Y1D) 

IY1D=AIMAG<Y1D)X 
RJ1D=REAL< J1D) 

I J1D=AIMAG ( J1D ) 

RY1D=REAL< Y1D) 

IY1D=AIMAG< Y1D) 

H1A=CMPLX(RJ1A+IY1A» I J1A-RY1A) 
H1D=CMPLX<RJ1D+IY1D»+I J1D-RY1D) 

GR1=-H1 A*  < -CMPLX <  BZ  t -AZ ) *AA/KR ) 

GR1 =GR1-H 1 D#  < -CMPLX  <  BZ » -AZ ) *AA*M/KR ) 
GR2=GR1#<0.0»-1.0) 
GRAl=GRl*CONJG<F)+CONJG<GRi)*F 
GRA2=GR2*C0N JG ( F ) +CON JG  <  GR2 ) *F 
GRAD ( 1 ) “REAL (GRAl) 

GRAD<2)=REAL(GRA2) 

RETURN 

END 
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Output  Saapl*  for  Prograa  No.  4 


D/A  -  10 

D  IS  t  .1000 

REAL  KZ  IS  J  -12.3694  IHA0  KZ  IS  I  68.1768 
THE  CURRENT  IS  J  0. 13376462E-03 
CUR/TEM  ISJ  0. 36961 400E-01 

D  IS  J  .2000 

REAL  KZ  IS  I  -6.2624  I MAG  KZ  IS  J  33.6657 

THE  CURRENT  IS  J  0.27066794E-03 
CUR/TEM  ISJ  0 . 74790076E-01 

D  IS  J  .3000 

REAL  KZ  IS  J  -4.2655  IMAG  KZ  IS  J  21.9671 
THE  CURRENT  IS  J  0.41419069E-03 
CUR/TEM  ISJ  0 . 1 1 444781 E+00 

D  IS  J  .4000 

REAL  KZ  IS  J  -3.3017  IMAG  KZ  IS  J  15.9635 
THE  CURRENT  IS  J  0 .56857456E--03 
CUR/TEM  ISJ  0. 15710665E+00 

»  IS  J  .  .'000 

REAL  KZ  IS  J  -2.758 7  IMAG  KZ  IS  t  12*2276 

THE  CURRENT  IS  •  0.73941571E-03 
CUR/TEM  ISJ  0 . 20431292E+00 

D  IS  J  ,6000 

REAL  KZ  IS  J  -2.4367  IMAG  KZ  IS  J  9.6135 

THE  CURRENT  IS  J  0.93456462E-03 
CUR/TEM  ISJ  0 . 25823583E+ 00 
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THE  PURP08E3  OF  THX8  PROGRAM  ARES 

1 >  TO  FIND  THE  REAL  AND  IHAOINARY  PARTS  OF  THE  R-DIRECTED 
AND  Z-DXRECTED  NAVE  NUMBERS » ASSOCIATED  HSTH  THE  HZOHER 
ORDER  EVEN  LEAKY  MODES  IN  THE  CA8E  OF  VOLT  EXCITATION 
2)  TO  CALCULATE  THE  CONTRIBUTION  OF  EACH  MODE  TO  THE  TOTAL 
CURRENT* 

THE  PROGRAM  VARIES  THE  RATXO-^D/A*  * WHILE  KEEPING 
THE  DISTANCE  BETWEEN  THE  WXRE8-*Da  C0N8TANT 


COMMON  KR'AA»H»DR;AR 
COMPLEX  KRrKZP'KRP'RES 

REAL  I JO A  t  X JODt  X J1A»  X J1D»  XYOAt  XYODt  XY1A» IYID 
COMPLEX  KRA  t KRD » JOA » YO A  » JOD » YOD » HOA  >  HOD  »  J1 A 

COMPLEX  Y1A»  JlD'YlD»GRl».QR2tHlA»HlDtFfGRAl»GRA2 
COMPLEX  KRPAtKRPD 
REAL  KZ(2>  »K*KH 
REAL  GRAD  <  2  > 

REAL  0GRAD<2)»C0N<2> 

DEPARTURE  POINT  FOR  THE  SEARCHING  XN  THE  KZ  PLANE 

1-1 

BR-1. 

AR— 2. 

START  SEARCHING  FOR  THE  R00T8  XN  THE  KZ  PLANE 

IF  (BR*GT » IOO)  GO  TO  8 
WRITE  <5'18>  I 

FORMAT <///» IX » 'THE  INDEX  OF  THE  MODE  IS  '»I2»//> 

THE  DISTANCE  BETWEEN  THE  WIRES 

D-2.0 

WRITE  <5»80)  D 

FORMAT <////» IX* '  D  IS  '*F6*4*///> 

INITIAL  VALUE  FOR  THE  RATIO  "D/A* 


M-10 

IF(M*GT « 1000)  GO  TO  60 
0GRAD(1 )“0. 

0Gi*AD<2>-0. 

COM  1  >-0.50 
C0N<2)-0*S0 
N-2 

PI-3* 1415927 
K-2SPI 

BZ-SQRT <  0 . 5* ( KH+SQRT ( ( KH**2+ ( 2SARSBR ) t*2 ) ) > ) 
AZ— ARSBR/BZ 
KR-CMPLX  <  BR  * -AR ) 

KZ(1)-BZ 
KZi2)-AZ 
AA-1 */M 
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noon 


WRITE  <5»19)  M 
FORMAT <//» IX r'B/A  •#*M> 

CALL  FUNCT<NtKZ»VAL.ORAD> 

XF<VAL»LT«0.lE-S>  00  TO  1 
IF  (QRAD(l))  22»2?»27 
IF  <06RAD<1>>  23»24»24 
CQNd  )*CON(l)/2« 

KZd>«KZ(l>+CONU> 

OGRAD d ) »GRAD d ) 

60  TO  31 

IF  (OORAD(l) )  28*29*29 

C0Nd>-C0Nd>/2 

KZd>«KZCl>-CON<l> 

06RAD ( 1 ) "GRAD ( 1 ) 

CONTINUE 

IF  <GRAD<2>>  32*37*37 
IF  <0GRAD<2) >  33*34*34 
CON<2)-CON<2>/2« 

KZ  <  2 )  *KZ  ( 2 )  4-CON  (  2  > 

OORAD ( 2  >  *GRAD ( 2  > 

GO  TO  3 

IF  <0GRAD<2>>  38*39*39 
CON ( 2  > "CON  <  2  > /2 
KZ ( 2  >  «KZ ( 2 ) -CON ( 2  > 

OGRAD <  2  > “GRAD <  2 ) 

GO  TO  3 

CONTINUE 

A-D/N 

BP-+BR/D 

AP-+AR/D 

KRP»CMPLX<BP*-AP> 

KRPA-KRP9A 

KRPD-KRP4D 

K-2*PI 

KH«K**2+AP**2-BP*t2 

BZ-SQRT  <  0 . 3*  C KH+8QRT <  <  KH*t2+ ( 28AP8BP >**2) >>> 

AZ—APBBP/BZ 

KZP»CMPLX(-BZ*AZ> 

WRITE (3, 90)  KZP 

FORMAT</ *1X* "REAL  KZ  18  l',F9.4*'  IMAG  KZ 
CALL  CBE81 <  KRPA  *  J1 A *  Y1 A  > 

CALL  CBE81 (KRPO  *  J1D » Y1D ) 

RJIA-REAL(JIA) 

I J1A*AIMAG( J1A) 

RY10«REAL(Y1D) 

IY1D-AIMAG(Y1D) 

RJID-REAL(JID) 

IJ1D-AXNAG(J1D> 

RYID-REAL(YID) 

lYlO-AXMAG(YlD) 

H1A-CMPLX(RJ1A+IY1A» IJ1A-RY1A) 

H1D-CMPLX  <  RJ1D+IY1D * +1 J1D-RY1D > 

THE  CONTRIBUTION  OF  THE  LEAKY  MODE 
TO  THE  TOTAL  CURRENT  NORMALIZED  TO  THE 
CONTRIBUTION  OF  THE  TEM  MODE 


IMAG  KZ  18  t  *  »F9*4) 
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c 

RES-4./377. 

RES*  RES* ( 0 1 1 > *K/KRP/KZP/ ( A*H1 A+DtHID > 

ARS*CABS<R£S> 

WRITE <3» 40 >  ARS 

40  F0RHAT</»1X»'TH£  CURRENT  18  t'»E15*B> 

T-M 

RAT*AR8*377 ♦ *AL0G  <  T  > /PI 
WRITE<5*47>  RAT 

47  F0RHAT</»1X» 'THE  CUR/TEH  IS  l',E15.B> 

C 

C  CHANGING  THE  RATIO  *D/A* 

C 

M»M+1 
GO  TO  41 
40  CONTINUE 

C 

C  CHANGING  THE  DEPARTURE  POINT  IN  ORDER  TO 

C  START  THE  SEARCHING  FOR  THE  CONSECUTIVE  ROOT 

C 

I-I+l 
BR-BR+5 
GO  TO  9 

8  CONTINUE 

STOP 
END 
C 

C  THE  SUBROUTINE  CBESO  CALCULATES  THE  ZERO  ORDER  BESSEL 

C  FUNCTIONS  FOR  A  GIVEN  INPUT  VARIABLE  Z 

C 

SUBROUTINE  CBESO(Z»BSJO»BSYO> 

COMPLEX  Z»BSJO»BSYO»YrW»PO>FOf CEXP 

COMPLEX  CCOS iCSINn CSGRT > CLOG r SPO r  SQO  »  QO  f  CS » SN , SP 1 » SC 1 

PI=3. 141593 

A=CABS<Z> 

BSJO=( 1 . »0 .  ) 

IF(A.EQ.O. )  GO  TO  2 
IF <  A.GT .3.  )  GO  TO  1 
Y=Z#Z/9. 

BS J0=1 . +Y* ( -2 ♦ 2499997+Y*  < 1 . 2656208+Y* ( - ♦3143846 
1 +Y* < .0444479+Y*  < - . 0039444+Y* . 00021 ) > ) ) ) 

GO  TO  2 

1  U*3./Z 

F0» . 79788456+W*  < - . 00000077+U* ( - . 0055274+W*  < - . 000095 1 2 
1 +W*  < . 00 137237+U* ( - . 00072805+W* ♦ 00014474 ) ) > ) > 

P0= ♦ 78539814+W*( ♦ 04166397+W*  C . 00003954+W*<- . 00262573 
1+W*< .00054125+U*< ♦00029333-U*. 00013558 >  > ) ) ) 

BS J0«F0*CC0S ( Z-PO ) /CSQRT <  Z ) 

2  CONTINUE 
BSYO=<-l.E30tO> 

IFCA.EQ.O, )  GO  TO  3 
IF(A,GT.3.)  GO  TO  4 

BSYO* ♦ 63661977*CL0G( ♦ 5#Z ) *BS J0+ ♦ 36746691+Y*  < .60559366 
l+Y*<-.74350384+Y*( .25300117+Y*(-.04261214+Y*< .00427916 
2-Y# ♦ 00024846 ) ) ) ) ) 

GO  TO  3 
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CONTINUE 

BSYO-FOBCS IN ( Z-PO ) /CSQRT ( Z  > 

3  CONTINUE 
RETURN 
END 

THE  SUBROUTINE  CBE81  CALCULATES  THE  FIR8T  ORDER  BE88EL 
FUNCTIONS  FOR  A  GIVEN  INPUT  VARIABLE  Z 

SUBROUTINE  CBEB1<Z»BSJ1»B8Y1> 

COMPLEX  Z'B8J1'BSY1*Y»U'P1'F1»CEXP 

COMPLEX  CCOS  t CSQRT » CLOG  t CSIN  t SP1 » 801 1 Gi *  F » CS » 8N » 8P2 » SQ2 
PI-3.141593 
A-CABS(Z) 

BSJ1-0. 

IF( A.EQ.O.  >  GO  TO  2 
IFCA.GT .3.  >  GO  TO  1 
Y-ZBZ/9. 

BSJ1-Z*< . 5+Y*  < - ♦ 562499B5+Y* ( « 21093573+YB ( - » 03954289 
1+Y*( » 00443319+Y*  < - .00031 761+YB . 00001109  > ) > ) > ) 

GO  TO  2 
U-3./Z 

Fl-.79788456+U*< .00000136+U*< . 01659A 67+U*< .00017105 
1+W*<- .0024951 1+W*< .00 113653-U*. 00020033) >>) ) 

PI -2 . 3561 9449+ W*  < - . 124996 12+W* ( - . 00005650+U* ( . 00437879 
1 +W*  < - . 00074348+W# ( - . 00079B24+W* . 00029166  >  >  >  > ) 

BS J1-F14CC0S <  Z-Pl ) /CSQRT ( Z  > 

CONTINUE 
BSY1*<-1 .E30»0) 

IF< A.EQ.O. >  GO  TO  3 
IF<A.GT,3.)  GO  TO  4 

BSYl-<-.63661977+Y*< . 2212091+Y*<2. 1682709+Y*< -1 .3164827 
1+Y*< .312395i+Y*<-. 0400976+Y*. 0027873) > ) ) ) )/Z 
2+ ♦ 636619774CL00 ( .5*Z)*BSJ1 
GO  TO  3 

BSY1»F1*CSIN<Z-P1)/CSGRT<Z) 

CONTINUE 
RETURN 
END 

THE  SUBROUTINE  FtINCT  CALCULATES  FOR  EACH  VECTOR  KZ 
THE  VALUE  OF  THE  FUNCTIONAL  WHICH  WE  TRY 
TO  MINIMIZE  AND  ITS  GRADIENT 

SUBROUTINE  FUNCT <  N » KZ » VAL » GRAD ) 

COMMON  KR'AA'MfBR'AR 
DIMENSION  GRAD<2> 

REAL  KZ<2) »K»KH 

REAL  I JOA  » I J0D» I J1A» IJ1D»IY0A»IY0D»IY1A»IY1D 
COMPLEX  KR . KR A » KRD » JOA  »  YOA  t JOD  » YOD  » HOA » HOD  >  J1 A 
COMPLEX  YlA»JlD?YlDrGRl *GR2»H1A»H1D»F*GRA1 »GRA2 
BZ-KZ< 1 ) 

AZ»KZ<2) 

PI-3.1413927 
K-2*PI 

KH«K**2+AZ**2-BZ**2 
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BR«8QRT ( 0 » St ( KH+8GRT < <KH*«2+<2*AZ*BZ>**2> ) >  > 

AR—AZBBZ/BR 

KR-CMPLX<BR»-AR> 

KRA-KRBAA 

KRD*KRA*H 

CALL  CBESO(KRA» JOAtYOA) 

CALL  CBE80(KRD> JODtYOD) 

RJOA-REAL(JOA) 

XJOA-AZHAO(JOA) 

RYOA-REAL(YOA) 

XYOA*AXHAG(YOA) 

R JOD-REAL ( JOB ) 

I JOD*AXMAG( JOD) 

RYOD-REAL(YOD) 

XYOD-AXMAG(YOD) 

HO A-CHPLX ( RJOA+ X  YOA  »  X  JOA-RYOA > 

HOD-CMPLX <  RJOD+ 1 YOD 1 1 JOD-R YOD > 

F-HOA+HOD 

VAL-CABS(F) 

VAL«VAL**2 

CALL  CBES1 (KRA» J1A» Y1A> 

CALL  CBES1 (KRDf JIDt Y1D) 

RJlA»REALt J1A) 

IJIA-AIMAG(JIA) 

RY1D«REAL<Y1D> 

IY1D*AIMAG< Y1D) 

RJ1D-REALCJ1D) 

IJ1D«AIHAG<J1D> 

RY1D»REAL<Y1D> 

IY1D«AIMAG<Y1D> 

HlA»CMPLX(RJiA+IYlAr I J1A-RY1A) 
H1D*CMPLX<RJ1D+IY1Df+I J1D-RY1D) 

OR  1 —HI  A*  <  -CMPLX  <  BZ » -AZ  >  BAA/KR  ) 

GR 1 »GR 1 -H 1 DB  < -CHPLX <  BZ  f - AZ  > * AABM/KR ) 
GR2»GR1*(0«0f-1.0) 
GRA1»GR1BC0NJG(F)+C0NJG(GR1 >BF 
GRA2=GR2BCQN JG ( F ) +CON JG ( GR2 >  BF 
GRAD ( 1 ) “REAL ( GRA1 ) 

GRAD ( 2 ) "REAL ( GRA2  > 

RETURN 

END 
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D/A  «  10 

REAL  KZ  IS  X  -6,0968  I HAG  KZ  IS  t  0.1S82 

THE  CURRENT  IS  X  0.23093203E-02 
THE  CUR/TEH  IS  X  0.63810379E+00 

D/A  »  11 

REAL  KZ  IS  :  -6,1045  IMAQ  KZ  IS  t  0,1595 

THE  CURRENT  IS  X  0 . 22558292E-02 
THE  CUR/TEM  IS  X  0.64912435E+00 

D/A  *  12 

REAL  KZ  IS  ‘,  -6,1111  IMAG  KZ  IS  X  0,1605 

THE  CURRENT  IS  X  0.22111599E-02 
THE  CUR/TEM  IS  X  0 . 65935864E+00 

D/A  =  13 

REAL  KZ  IS  J  -6,1169  I HAG  KZ  IS  :  0.1616 

THE  CURRENT  IS  X  0.21676725E-02 
THE  CUR/TEM  IS  :  0.66721214E+00 

D/A  *  14 

REAL  KZ  IS  X  -6.1220  IMAG  KZ  IS  X  0,1626 

THE  CURRENT  IS  X  0.21300598E-02 
THE  CUR/TEM  IS  :  0.67457789E+00 

D/A  *  15 


REAL  KZ  IS  :  6.1265  IMAG  KZ  IS  X  0.1636 


Program  No.  6 


C  .  . . . . . . . 

C  THE  PURPOSES  OF  THIS  PROGRAM  ARE* 

C  1)  TO  FIND  THE  REAL  AND  IMAGINARY  PARTS  OF  THE  R-DIRECTED 

C  AND  Z-DIRECTED  WAVE  NUMBERS f ASSOCIATED  WITH  THE  HIGHER 

C  ORDER  ODD  LEAKY  MODES  IN  THE  CASE  OF  VOLT  EXCITATION 

C  2)  TO  CALCULATE  THE  CONTRIBUTION  OF  EACH  MODE  TO  THE  TOTAL 

C  CURRENT. 

C  THE  PROGRAM  VARIES  THE  RATIO-'D/A* » WHILE  KEEPING 

C  THE  DISTANCE  BETWEEN  THE  WIRES-*D*  CONSTANT 

C  * . . . . 

COMMON  KR'AAfMrBR’AR 
COMPLEX  KR»KZP»KRP»RES 

REAL  IJCA  » I JODf IJlArl JIDflYOAf IY0D»IY1A» IY1D 
COMPLEX  KRA » KRD » JOA » YOA» JOD  » YOD  t HOA  »  HOD  »  J1 A 

COMPLEX  Y1A» JlD»YlD»GRi»GR2»HlAfHlD»F»GRAl»GRA2 
COMPLEX  KRPA»KRPD 
REAL  KZ(2) f Kf KH 
REAL  GRAD  <  2 ) 

REAL  OGR AD  <  2 ) » CON  <  2 ) 

C 

C  DEPARTURE  POINT  FOR  THE  SEARCHING  IN  THE  KZ  PLANE 

C 

1  =  1 
BR=1 . 

AR=-2. 

C 

C  START  SEARCHING  FOR  THE  ROOTS  IN  THE  KZ  PLANE 

C 

9  IF  (BR.GT.lOO)  GO  TO  8 

WRITE  <5»lu>  I 

18  FORMAT (///»1X»' THE  INDEX  OF  THE  MODE  IS  'rI2r//> 

C 

C  THE  DISTANCE  BETWEEN  THE  WIRES 

C 

D=2 . 0 

WRITE  <5f80)  D 

SO  F0RMAT(////»1X» '  D  IS  '»F6.4f///> 

C 

C  INITIAL  VALUE  FOR  THE  RATIO  "D/A* 

C 

M=10 

41  IF <M.GT. 1000)  GO  TO  60 

OGRAD< 1 >=0 . 

0GRAD<2>=0. 

CON< 1 ) =0.50 
CCN(2>=0 . 50 
N=2 

PI=3. 1415927 
K»2*PI 

KH»K**2+AR**2-BR**2 

B^»SORT  ^  0 . 5*  <  KH+SQRT ( <  K H**2+ <  2*AR*BR ) **2  >  >  > ) 

AZ»-AR#BK/BZ 
KR=CMPLX  <  BR » -AR  > 

KZ( 1 )=BZ 
KZ(2)*AZ 
AA-l./M 


WRITE  <5*19>  N 

19  FORMAT (// »1X» 'D/A  -'*I4> 

3  CALL  FUNCT<N*KZ*VAL»0*AD) 

IF(VAL.LT.0.1E-5>  00  TO  1 
ZF  (ORADd  >  >  22*27*27 

22  IF  (OORAD(l) )  23*24*24 

24  CON ( 1 ) *CON ( 1 > /2 » 

23  KZ ( 1 > “KZ d  >  +CON ( 1 > 

OGRAD(l)  -ORADd) 

00  TO  31 

27  IF  ( 06RAD < 1 ) >  28*29*29 

28  CON  < 1 ) “CON d  > /2 

29  KZ < 1 ) “KZ d  > -CON Cl) 

OGRAD  (1 )  “GRAD  (1  > 

31  CONTINUE 

IF  (GRAD<2) )  32*37*37 

32  IF  (C0RAD<2) )  33*34*34 

34  CON ( 2 ) =CON ( ?  > /2 « 

33  KZ ( 2 ) =KZ ( 2 ) *CON <  2 ) 

OGRAD ( 2  >  =GRAD ( 2 ) 

GO  TO  3 

37  IF  <0GRAD<2>>  38*39*39 

38  CON  <  2 ) =CON  <  2 ) /2 

39  KZ ( 2 ) “KZ ( 2 ) -CON  <  2  > 

OGRAD ( 2 ) =GRAD ( 2 ) 

GO  TO  3 

1  CONTINUE 

70  A“D/M 

BP=+BR/D 

AP=+AR/D 

KRP“CMPLX (BP  * -AP ) 

KRPA“KRP*A 

KRPD=KRP*D 

K“2*PI 

KH=K**2+AP**2-BP**2 

BZ=SQRT<0.5*(KH+SQRT(  (i>H**2+(2*AP*BP)**2>  > ) ) 

AZ=-AP*BP/BZ 

KZP=CMPLX<-BZ»AZ) 

WRITE(5*90>  KZP 

90  FORMAT (/ *  IX* 'REAL  KZ  IS  J,*F9.4*/  IMAG  KZ  IS  J'*F9.4) 

CALL  CBES1 (KRPA* J1A*Y1A) 

CALL  CBES1 <KRPD» J1D* Y1D) 

RJ1A=REAL ( J1A  > 

IJ1A=AIMAG< J1A> 

RY1D“REAL<Y1D) 

IY1D=AIMAG<Y1D> 

RJ1D“REAL ( J1D ) 

I J1D=AIMAG( J1D) 

RY1D“REAL  <  Y1D ) 

IY1D“AIMAG( Y1D) 

H1A=CMPLX<RJ1A+IY1A* I J1A-RY1A) 

H1D=CMPLX  (RJ1D+I Y1D*  +  1  J1D-RY1.D) 

C 

C  THE  CONTRIBUTION  OF  THE  LEAKY  MODE 

C  TO  THE  TOTAL  CURRENT  NORMALIZED  TO  THE 

C  CONTRIPUT’IN  OF  THE  TEM  MODE 


c 

RES-4./377. 

RES*  RES*(Orl)*K/KRP/KZP/(A*HlA+D*HlD> 

ARS-CABS(RES) 

WRITE (3, 40)  ARS 

40  FORMAT (/» IX » 'THE  CURRENT  ZS  t'rE15.8) 

T*H 

RAT*ARS*377 . *ALOG  <  T  > /PI 
WRITE<5»47)  RAT 

47  FORMAT </» lXt 'THE  CUR/TEM  IS  :'rE15.8> 

C 

C  CHANGING  THE  RATIO  *D/A* 

C 

M*M+1 
GO  TO  41 
60  CONTINUE 

C 

C  CHANGING  THE  DEPARTURE  POINT  IN  ORDER  TO 

C  START  THE  SEARCHING  FOR  THE  CONSECUTIVE  ROOT 

C 

1=1+1 
BR=BR+5 
GO  TO  9 

8  CONTINUE 

STOP 
END 
C 

C  THE  SUBROUTINE  CBESO  CALCULATES  THE  ZERO  ORDER  BESSEL 

C  FUNCTIONS  FOR  A  GIVEN  INPUT  VARIABLE  Z 

C 

SUBROUTINE  CBESO < ZrBS JO rBSYO) 

COMPLEX  Z * BS JO fBSY0»Y»W»P0» FOr CEXP 

COMPLEX  CCOS  r  CS I N » CSQRT » CLOG  r SPO  r  SQO » QO  »  CS » SN » SP1 » SQ 1 

PI=3. 141593 

A=CABS(Z) 

BSJ0=<1. .0. ) 

IF(A.EQ<0. )  GO  TO  2 
IF<A»GT.3. )  GO  TO  1 
Y=Z*Z/9. 

BS  J0*1 .  +Y*  ( -2 . 2499997+Y*  <  1 . 2656208+Y*  <  -  .3163866 
1+Y*< .0444479+ Y*(-.0039444+Y*. 00021 >  >  > )  > 

GO  TO  2 

1  W=3./Z 

F0= . 797BB456+W*  < - . 00000077+W* ( - .0055274+W*  < - . 00009512 
1+W*< .00137237+W*<-. 00072805+W*. 00014476) > > ) ) 

P0= ♦ 7853981 6+W#  < ,04166397+W*C ,00003954+W*<-. 00262573 
1+W#< .00054125+W*< .00029333-W*. 00013558) ) >  >  > 

BS J0=F0*CC0S <  Z-PO  > /CSQRT  <  Z  > 

2  CONTINUE 
BSY0=<-1.E30>0) 

IF (A.EQ.O. )  GO  TO  3 
IF< A.GT.3. )  GO  TO  4 

BSY0= ♦ 63661977*CL0G ( . 5*Z)*BSJ0+ .36746691 +Y# ( . 60559366 
l+Y#<-.74350384+Y#<  «25300117+Y$(-.04261214+Y*( .00427916 
2- Y*. 00024846) ) ) >) 

GO  TO  3 
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CONTINUE 

BS  YO*FO*CSIN  <  Z-PO  >  /CSQRT  <  Z  > 

3  CONTINUE 
RETURN 
END 

C 

C  THE  SUBROUTINE  CBES1  CALCULATES  THE  FIRST  ORDER  BES8EL 

C  FUNCTIONS  FOR  A  GIVEN  INPUT  VARIABLE  Z 

C 

SUBROUTINE  CBESl<ZtB8Jl»BSYl> 

COMPLEX  ZrBSJl»BSYl»Y»W»Pl,Fl»CEXP 

COHPLE X  CCOS  *  CSQRT  »  CLOG  t  CSIN  »  SP1 »  SOI r 01 »  F  t CS  t SN  r  SP2  » SQ2 

PI-3. 141593 

A«CABS(Z> 

BSJl-O. 

IF < A.EQ.O. >  GO  TO  2 
IF( A.GT.3* )  GO  TO  1 
Y*Z*Z/9. 

BSJ1*Z*< •  5+ Y#  < - . 56249985+Y*  < .21093573+Y*<-. 03954289 
1+Y*< ,00443319+Y*<-. 00031761+Y*. 00001109) > ) > ) ) 

GO  TO  2 

1  W=3./Z 

F 1 = ♦ 79788456+W* ( . 00000156+U* ( ,01659667+W*  ( .00017105 
1+W*<-. 0024951 1+W*< . 00 1136S3-U*. 00020033 >  > >  > ) 

Pl=2 . 35619449+W* < - . 12499612+W*< - . 00005650+W* ( . 00637879 
l+W*(-.0007434B+U*<-.00079824+W*.00029166) >  > ) ) 

BS J1 =F 1 *CCOS ( Z-Pl ) /CSQRT ( Z  > 

2  CONTINUE 
BSY1=<-1 . E30.0) 

IF (A.EQ.O. )  GO  TO  3 
IF (A.GT .3. )  GO  TO  4 

BSY1=  < - . 63661977+Y*  < . 2212091+Y*  <  2 . 1682709+Y* ( -1 .3164827 
1+Y*<  ♦  3 123951 +Y* < -  •  0400976+Y* . 0027873 )  )  ) )  >/Z 
2+.63661977*CL0G( .5*Z)*BSJ1 
GO  TO  3 

4  BSY1=F1*CSIN<Z-P1)/CSQRT(Z> 

3  CONTINUE 
RETURN 
END 

C 

C  THE  SUBROUTINE  FUNCT  CALCULATES  FOR  EACH  VECTOR  KZ 

C  THE  VALUE  OF  THE  FUNCTIONAL  WHICH  WE  TRY 

C  TO  MINIMIZE  AND  ITS  GRADIENT 

C 

SUBROUTINE  FUNCT ( N  »  KZ » VAL  *  GRAD ) 

COMMON  KR'AA'MfBR'AR 
DIMENSION  GRAD ( 2 ) 

REAL  KZ<2) »K»KH 

REAL  I JOA  »IJOD»IJlA»IJlD»IYOArIYOD>IYlA*IYlD 
COMPLEX  KR.KRA.KRDu JOAr YOAf JODr YODfHOAf HODr J1A 
COMPLEX  Y1A»J1DfY1D»GR1 »GR2»HlAf H1D»F »GRA1 >GRA2 
BZ=<KZ<1> 

AZ*KZ(2> 

PI®3* 1415927 
K=2*PI 

KH=K**2+AZ**2-BZ*#2 
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BR*3QRT <  0 . 3t  < KH+SGRT < <KH**2+<2*AZ*BZ>**2> >  > ) 

AR»-AZ*BZ/BR 

KR"=CMPLX  <  BR » -AR  > 

KRA-KRBAA 

KRD*KRA«M 

CALL  CBESO(KRA» JOA» YOA) 

CALL  CBESO(KRD»JOD»YOD> 

R JOA*REAL ( JOA > 

X JOA* A I MAG (JOA) 

RYOA*REAL  <  YOA  > 

IYOA»AIMAG<YOA) 

R JODaREAL ( JOD  > 

I JOD=AIMAG( JOD) 

RYODaREAL  <  YOD ) 

IYODaAIHAG(YOD> 

HOAaCMPLX ( RJOA+IYOA » I JOA-RYOA > 

HOD*CMPLX <  R JOD+I YOD  » I JOD-RYOD  > 

FaHOA-HOD 

VAL=CABS<F> 

VAL=VAL**2 

CALL  CBES1 <KRA» J1A>Y1A) 

CALL  CBES1 (KRD» J1D»Y1D) 

RJ1A«REAL(J1A> 

IJ1A*AIMAG<J1A> 

RYlDaREAL(YtD> 

IY1D=»AIMAG(  Y1D) 

RJ1D=REAL  <  J1D  > 

I J1D=AIMAG< J1D) 

RYlDaREAL(YlD) 

IY1D=AIMAG<Y1D> 

H1A=CMPLX (RJIA+IYlAr IJ1A-RY1A) 

HlDaCMPLX (RJID+IYIDr  +1 J1D-RY1D ) 

GR1 »-Hl A*  < -CMPLX <  BZ  » -AZ ) BAA/KR > 

GR 1  aGR  1  -H 1 DB  < -CMPLX  C  BZ » - AZ  > B AABM/KR ) 
GR2=GRl»(O.Or-1.0> 
GRA1=GR1BCONJG(F)+CONJG<GR1 >*F 
GRA2=GR2*C0NJG <F  >  +CONJG < GR2 ) *F 
GRAD< l>aREAL <GRA1> 

GRAD ( 2 ) =REAL  <  GRA2 ) 

RETURN 

END 
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Output  Staple  for  Program  Ho.  6 


D  IS  2*0000 


D/A  ~  10 

REAL  KZ  IS  J  -5.3197  I MAG  KZ  IS  : 
THE  CURRENT  IS  i  0 . 17374900E- 02 
THE  CUR/TEM  IS  X  0 . 48009754E+00 

D/A  =  11 

REAL  KZ  IS  i  -5.3473  IMAG  KZ  IS  X 
THE  CURRENT  IS  X  0 . 16872877E-02 
THE  CUR/TEM  IS  X  0 . 48552414E+00 

D/A  «  12 

REAL  KZ  IS  X  -5.3704  IMAG  KZ  IS  { 
THE  CURRENT  IS  X  0 . 16398675E-02 
THE  CUR/TEM  IS  {  0. 489001 64E+00 

D/A  =  13 

REAL  KZ  IS  X  5.3899  IMAG  KZ  IS  } 
THE  CURRENT  IS  X  0 . 15975309E-02 
THE  CUR/TEM  IS  }  0.49172189E+00 

D/A  =  14 

REAL  KZ  IS  {  -5.4066  IMAG  KZ  IS  X 

THE  CURRENT  IS  }  0 . 15589255E-02 
THE  CUR/TEM  IS  X  0 . 49370289E +00 

D/A  =  15 

REAL  KZ  IS  X  -5,4212  IMAG  KZ  IS  X 
THE  CURRENT  IS  {  0 . 15255773E-02 
THE  CUR/TEM  IS  X  0.49577245E+00 


0.3963 


0.3998 


0.4036 


0.4071 


0.4106 


0.4135 
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